Plotting 1 arc-minute global relief map in Python (codes included)

Utpal Kumar   1 minute read      

Plotting 1 arc-minute topographic map in python

  • see here too
  • uses Python’s matplotlib (basemap) library for plotting
  • to use topographic map from 1 arc-minute global relief data from NOAA (or even higher resolution)
  • uses a sliced array for plotting the data for the selected region
  • run-time has been optimized
  • colormap is customizable
  • visit here for the codes
  • 1 arc-minute global relief data can be downloaded from NOAA’s webpage.

World Map

## World map
lonmin, lonmax = -180,180
latmin, latmax = -50,50

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
map = Basemap(projection='merc',resolution = 'l', area_thresh = 10000., llcrnrlon=lonmin, llcrnrlat=latmin,urcrnrlon=lonmax, urcrnrlat=latmax)
cs = plot_topo(map,cmap='terrain',lonextent=(lonmin, lonmax),latextent=(latmin, latmax),zorder=2)
fig.colorbar(cs, ax=ax, shrink=0.4)
map.drawcoastlines(color='k',linewidth=0.5)
map.drawcountries(color='k',linewidth=0.1)
# map.drawstates(color='gray',linewidth=0.05)
# map.drawrivers(color='blue',linewidth=0.05)
parallelmin = int(latmin)
parallelmax = int(latmax)+1
map.drawparallels(np.arange(parallelmin, parallelmax,10,dtype='int16').tolist(),labels=[1,0,0,0],linewidth=0,fontsize=6)

meridianmin = int(lonmin)
meridianmax = int(lonmax)+1
map.drawmeridians(np.arange(meridianmin, meridianmax,20,dtype='int16').tolist(),labels=[0,0,0,1],linewidth=0,fontsize=6)

map.drawmapboundary(color='k', linewidth=2, zorder=1)

plt.savefig('world_map.png',bbox_inches='tight',dpi=300)
plt.close('all')

world map

Africa Map

## Map1
lonmin, lonmax = -30,70
latmin, latmax = -40,40

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
map = Basemap(projection='merc',resolution = 'l', area_thresh = 10000., llcrnrlon=lonmin, llcrnrlat=latmin,urcrnrlon=lonmax, urcrnrlat=latmax)
cs = plot_topo(map,cmap='terrain',lonextent=(lonmin, lonmax),latextent=(latmin, latmax),zorder=2)
fig.colorbar(cs, ax=ax, shrink=0.9)
map.drawcoastlines(color='k',linewidth=0.5)
map.drawcountries(color='k',linewidth=0.1)
map.drawstates(color='gray',linewidth=0.05)
map.drawrivers(color='blue',linewidth=0.05)
parallelmin = int(latmin)
parallelmax = int(latmax)+1
map.drawparallels(np.arange(parallelmin, parallelmax,10,dtype='int16').tolist(),labels=[1,0,0,0],linewidth=0)

meridianmin = int(lonmin)
meridianmax = int(lonmax)+1
map.drawmeridians(np.arange(meridianmin, meridianmax,20,dtype='int16').tolist(),labels=[0,0,0,1],linewidth=0)

map.drawmapboundary(color='k', linewidth=2, zorder=1)

plt.savefig('africa.png',bbox_inches='tight',dpi=300)
plt.close('all')
africa map

Africa Etopo (NOAA):

lonmin, lonmax = -30,70
latmin, latmax = -40,40

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
map = Basemap(projection='merc',resolution = 'l', area_thresh = 10000., llcrnrlon=lonmin, llcrnrlat=latmin,urcrnrlon=lonmax, urcrnrlat=latmax)
cs = plot_topo_netcdf(map,etopo_file='ETOPO1_Bed_g_gmt4.grd',cmap='terrain',lonextent=(lonmin, lonmax),latextent=(latmin, latmax),zorder=2)
fig.colorbar(cs, ax=ax, shrink=0.4)
map.drawcoastlines(color='k',linewidth=0.5)
map.drawcountries(color='k',linewidth=0.1)

parallelmin,parallelmax = int(latmin), int(latmax)+1
map.drawparallels(np.arange(parallelmin, parallelmax+1,10,dtype='int16').tolist(),labels=[1,0,0,0],linewidth=0,fontsize=6)
meridianmin,meridianmax = int(lonmin),int(lonmax)+1
map.drawmeridians(np.arange(meridianmin, meridianmax+1,20,dtype='int16').tolist(),labels=[0,0,0,1],linewidth=0,fontsize=6)
map.drawmapboundary(color='k', linewidth=2, zorder=1)
plt.savefig('test_africa.png',bbox_inches='tight',dpi=600)
plt.close('all')
africa map

Taiwan Plot

### ETOPO TAIWAN
lonmin, lonmax = 119,123
latmin, latmax = 20,26

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
map = Basemap(projection='merc',resolution = 'f', area_thresh = 1000., llcrnrlon=lonmin, llcrnrlat=latmin,urcrnrlon=lonmax, urcrnrlat=latmax)
cs = plot_topo_netcdf(map,etopo_file='ETOPO1_Bed_g_gmt4.grd',cmap='terrain',lonextent=(lonmin, lonmax),latextent=(latmin, latmax),zorder=2)
fig.colorbar(cs, ax=ax, shrink=0.9)

map.drawcoastlines(color='k',linewidth=0.5)
map.drawcountries(color='k',linewidth=0.1)
parallelmin,parallelmax = int(latmin), int(latmax)+1
map.drawparallels(np.arange(parallelmin, parallelmax+1,10,dtype='int16').tolist(),labels=[1,0,0,0],linewidth=0,fontsize=6)
meridianmin,meridianmax = int(lonmin),int(lonmax)+1
map.drawmeridians(np.arange(meridianmin, meridianmax+1,20,dtype='int16').tolist(),labels=[0,0,0,1],linewidth=0,fontsize=6)

map.drawmapboundary(color='k', linewidth=2, zorder=1)
plt.savefig('taiwan_plot.png',bbox_inches='tight',dpi=600)
plt.close('all')
africa map

See more: iescoders post

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.


Leave a comment