Topographic grid generation from etopo1

In this notebook, we will see how to generate a regular topographic grid from a particular region that will have the desired format for running a badlands simulation.

Badlands imports a regular grid that is later triangularized and used to compute erosion and deposition induced by a combination of climate, tectonics and sea-level changes.

Here we use ETOPO1 to create this initial grid. ETOPO1 is a 1 arc-minute global relief model of Earth's surface that integrates land topography and ocean bathymetry. It was built from numerous global and regional data sets, and is available in "Ice Surface" (top of Antarctic and Greenland ice sheets) and "Bedrock" (base of the ice sheets) versions. The dataset is made available from NOAA webserver.



In [1]:

    
%matplotlib inline

# Import badlands grid generation toolbox
import pybadlands_companion.toolGrid as tools

# display plots in SVG format
%config InlineBackend.figure_format = 'svg'

1. Define region bounding box

We first need to provide the extent of the area we want to get the topography from.

This is done by specifying the min/max longitudes and latitudes:

llcrnrlon: lower left corner longitude.
llcrnrlat: lower left corner latitude.
urcrnrlon: upper right corner longitude.
urcrnrlat: upper right corner latitude.

For more information regarding the function uncomment the following line.



In [2]:

    
#help(tools.toolGrid.__init__)



In [3]:

    
grid = tools.toolGrid(llcrnrlon = 173.5, llcrnrlat = -42.5, urcrnrlon = 174.5, urcrnrlat = -41.5)

2. Display region in Basemap

Before extracting the topography for the desired region, we can first add an ARCGIS image of the area to Basemap. We retreive the image directly from ArcGIS.com web map servers.

You can either use the EPSG registry to get a specific code for your region of interest or you can use the default EPSG code provided in this notebook (EPSG::3857).

EPSG::3857 is a projected CRS WGS 84/Pseudo-Mercator dataset which uses spherical development of ellipsoidal coordinates. The global data covers the world between 85.06°S and 85.06°N.

To visualise the Map using the EPSG code we use the plotEPSG function. For more information regarding the function uncomment the following cell.



In [7]:

    
#elp(grid.plotEPSG)



In [11]:

    
grid.plotEPSG( epsg=3857, llspace=0.25, fsize=(8,8), title = 'Map EPSG::3857 NE region Southern Island NZ' )

3. Extract dataset from etopo1

Instead of downloading the dataset on our local container, we use the NOAA THREDDS protocol to load the data from their server. The dataset is available in different format, here we choose to load the NETCDF file. Once the longitude and latitude for the entire grid has been loaded we only extract a subset (getSubset function) to download the values corresponding to the region of interest.

The getSubset function takes 3 parameters:

tfile: type of topographic / bathymetric file to load using THREDDS protocol (only etopo1 for now).
offset: offset to add to the grid to ensure the entire region of interest is still within the simulation area after reprojection in UTM coordinates.
smooth: use or not of a smoothing filter.

For more information regarding the function uncomment the following cell.



In [ ]:

    
#help(grid.getSubset)



In [ ]:

    
grid.getSubset(tfile = 'etopo1', offset = 0.1, smooth = False)

We know can use the mapUTM function to transform the dataset in UTM coordinates which is used to create Badlands grid.

In addition the function creates a plot of the region based on exported ETOPO1 dataset.

For more information regarding the function uncomment the following cell.



In [ ]:

    
#help(grid.mapUTM)



In [ ]:

    
grid.mapUTM(contour=50, fsize=(8,8), saveFig=False, nameFig='map')

4. Export regular grid data

In this final step, we export a CSV regular grid that will be loaded in the simulation. The file provides for each line the following information:

X coordinates in meters (this axis has a West to East orientation),
Y coordinates in meters (this axis has a South to West orientation),
Z coordinates in meters.

Note: Nodes must be defined in increasing order from the South/West corner, first along the X axis.

The buildGrid function takes the grid resolution from the ETOPO1 dataset and interpolates it on a new grid which has a resolution specified by the resolution parameter (in metres). The function provides 3 types of interpolation methods ('linear', 'nearest', 'cubic').

For more information regarding the function uncomment the following cell.



In [ ]:

    
#help(grid.buildGrid)



In [ ]:

    
grid.buildGrid(resolution=250., method='cubic', nameCSV='xyz1')

Using plotly, we visualise the created regular grid in 3D.



In [ ]:

    
#help(grid.viewGrid)



In [ ]:

    
grid.viewGrid(width=900, height=600, zmin=-5000, zmax=5000, title='Export Grid')



In [ ]: