[Python Ocean Topic 10] Cartopy draws terrain contour maps of specific areas
Ocean and Atmospheric Sciences
The previous issues can be considered to be about the drawing methods of plane elements
This issue is about plane painting in specific areas
Global regional bathymetry map
Contents of this issue
Drawing a floor plan of a specific area of an element: I have two methods:
The first one: Crop nc file
Second: Mask other areas and show only a specific area of the floor plan.
1: Crop nc file
# scs's range is lon from 100 to 123;lat from 0 to 25; print(np.where(lon >= 100))# 8400 print(np.where(lon >= 123))# 9090 print(np.where(lat >= 0))# 2700 print(np.where(lat >= 25))# 3450 lon1=lon[8400:9090] lat1=lat[2700:3450] ele1=ele[2700:3450,8400:9090]
Purpose Find the initial and ending positions of the range.
Large area proof data clipping!
2 Draw global areas but only show specific areas
ax.set_extent([100, 123, 0, 25], crs=ccrs.PlateCarree())# Set the display range
References and their role in this article
1: Python study notes: numpy selects qualified data: select, where, choose, nonzero – Hider1214 – Blog Park (cnblogs.com)
Full text code
1Global regional bathymetry map
# -*- coding: utf-8 -*-
# %%
# Importing related function packages
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as feature
import numpy as np
import matplotlib.ticker as ticker
from cartopy import mpl
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.font_manager import FontProperties
from netCDF4 import Dataset
from palettable.cmocean.diverging import Delta_4
from palettable.colorbrewer.sequential import GnBu_9
from paletable.colorbrewer.sequential import Blues_9
from paletable.scientific.diverging import Roma_20
from pylab import *
def reverse_colourmap(cmap, name='my_cmap_r'):
reverse = []
k = []
for key in cmap._segmentdata:
k.append(key)
channel = cmap._segmentdata[key]
data = []
for t in channel:
data.append((1 - t[0], t[2], t[1]))
reverse.append(sorted(data))
LinearL = dict(zip(k, reverse))
my_cmap_r = mpl.colors.LinearSegmentedColormap(name, LinearL)
return my_cmap_r
cmap = Blues_9.mpl_colormap
cmap_r = reverse_colourmap(cmap)
cmap1 = GnBu_9.mpl_colormap
cmap_r1 = reverse_colourmap(cmap1)
cmap2 = Roma_20.mpl_colormap
cmap_r2 = reverse_colourmap(cmap2)
# read data
a = Dataset('D:\pycharm_work\data\etopo2.nc')
print(a)
lon = a.variables['lon'][:]
lat = a.variables['lat'][:]
ele = a.variables['topo'][:,:]
lon1=lon[1:10800:110]
lat1=lat[1:5400:110]
ele1=ele[1:5400:110,1:10800:110]
print(len(lon1))
print(len(lat))
# Picture 3
#Set map global properties
scale = '50m'
plt.rcParams['font.sans-serif'] = ['Times New Roman'] # Set the overall font to Times New Roman
fig = plt.figure(dpi=300, figsize=(3, 2), facecolor='w', edgecolor='blue')#Set a drawing board and return it to fig
ax = fig.add_axes([0.05, 0.08, 0.92, 0.8], projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([0, 360, -90, 90], crs=ccrs.PlateCarree())# Set the display range
land = feature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
facecolor=feature.COLORS['land'])
ax.add_feature(land, facecolor='0.6')
ax.add_feature(feature.COASTLINE.with_scale('50m'), lw=0.3)#Add coastline: keyword lw sets line width; linestyle sets line style
cs = ax.contourf(lon1, lat1, ele1, levels=np.arange(-9000,0,20),extend='both',cmap=cmap_r1, transform=ccrs.PlateCarree())
# ------colorbar settings
cb = plt.colorbar(cs, ax=ax, extend='both', orientation='vertical',ticks=np.linspace(-9000, 0, 10))
cb.set_label('depth', fontsize=4, color='k')#Set the colorbar label font and size
cb.ax.tick_params(labelsize=4, direction='in') #Set the colorbar tick font size.
cf = ax.contour(lon, lat, ele[:, :], levels=np.linspace(-9000, 0, 5), colors='gray', linestyles='-',
linewidths=0.2,transform=ccrs.PlateCarree())
#Add title
ax.set_title('Etopo', fontsize=4)
# Use Formatter to format tick labels
ax.set_xticks(np.arange(0, 360, 30), crs=ccrs.PlateCarree())#Add latitude and longitude
ax.set_xticklabels(np.arange(0, 360, 30), fontsize=4)
ax.set_yticks(np.arange(-90, 90, 20), crs=ccrs.PlateCarree())
ax.set_yticklabels(np.arange(-90, 90, 20), fontsize=4)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.tick_params(color='k', direction='in')#Change the tick pointing inwards and set the color to blue
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, xlocs=np.arange(-180, 180, 30), ylocs=np.arange(-90, 90, 20),
linewidth=0.25, linestyle='--', color='k', alpha=0.8)#Add grid lines
gl.top_labels, gl.bottom_labels, gl.right_labels, gl.left_labels = False, False, False, False
plt.savefig('scs_elevation03.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1) # Output the map and set the border margin to be tight
plt.show()
2: Crop nc file;
# -*- coding: utf-8 -*-
# %%
# Importing related function packages
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as feature
import numpy as np
import matplotlib.ticker as ticker
from cartopy import mpl
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.font_manager import FontProperties
from netCDF4 import Dataset
from palettable.cmocean.diverging import Delta_4
from palettable.colorbrewer.sequential import GnBu_9
from paletable.colorbrewer.sequential import Blues_9
from paletable.scientific.diverging import Roma_20
from pylab import *
def reverse_colourmap(cmap, name='my_cmap_r'):
reverse = []
k = []
for key in cmap._segmentdata:
k.append(key)
channel = cmap._segmentdata[key]
data = []
for t in channel:
data.append((1 - t[0], t[2], t[1]))
reverse.append(sorted(data))
LinearL = dict(zip(k, reverse))
my_cmap_r = mpl.colors.LinearSegmentedColormap(name, LinearL)
return my_cmap_r
cmap = Blues_9.mpl_colormap
cmap_r = reverse_colourmap(cmap)
cmap1 = GnBu_9.mpl_colormap
cmap_r1 = reverse_colourmap(cmap1)
cmap2 = Roma_20.mpl_colormap
cmap_r2 = reverse_colourmap(cmap2)
# read data
a = Dataset('D:\pycharm_work\data\etopo2.nc')
print(a)
lon = a.variables['lon'][:]
lat = a.variables['lat'][:]
ele = a.variables['topo'][:]
print(lon)
print(lat)
# scs's range is lon from 100 to 123;lat from 0 to 25;
# so
print(np.where(lon >= 100))# 8400
print(np.where(lon >= 123))# 9090
print(np.where(lat >= 0))# 2700
print(np.where(lat >= 25))# 3450
lon1=lon[8400:9090]
lat1=lat[2700:3450]
ele1=ele[2700:3450,8400:9090]
# Picture 3
#Set map global properties
scale = '50m'
plt.rcParams['font.sans-serif'] = ['Times New Roman'] # Set the overall font to Times New Roman
fig = plt.figure(dpi=300, figsize=(3, 2), facecolor='w', edgecolor='blue')#Set a drawing board and return it to fig
ax = fig.add_axes([0.05, 0.08, 0.92, 0.8], projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([100, 123, 0, 25], crs=ccrs.PlateCarree())# Set the display range
land = feature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
facecolor=feature.COLORS['land'])
ax.add_feature(land, facecolor='0.6')
ax.add_feature(feature.COASTLINE.with_scale('50m'), lw=0.3)#Add coastline: keyword lw sets line width; linestyle sets line style
cs = ax.contourf(lon1, lat1, ele1, levels=np.arange(-9000,0,20),extend='both',cmap=cmap_r1, transform=ccrs.PlateCarree())
# ------colorbar settings
cb = plt.colorbar(cs, ax=ax, extend='both', orientation='vertical',ticks=np.linspace(-9000, 0, 10))
cb.set_label('depth', fontsize=4, color='k')#Set the colorbar label font and size
cb.ax.tick_params(labelsize=4, direction='in') #Set the colorbar tick font size.
cf = ax.contour(lon, lat, ele[:, :], levels=[-5000,-2000,-500,-300,-100,-50,-10], colors='gray', linestyles= '-',
linewidths=0.2,transform=ccrs.PlateCarree())
#Add title
ax.set_title('Etopo', fontsize=4)
# Use Formatter to format tick labels
ax.set_xticks(np.arange(100, 123, 4), crs=ccrs.PlateCarree())#Add latitude and longitude
ax.set_xticklabels(np.arange(100, 123, 4), fontsize=4)
ax.set_yticks(np.arange(0, 25, 2), crs=ccrs.PlateCarree())
ax.set_yticklabels(np.arange(0, 25, 2), fontsize=4)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.tick_params(color='k', direction='in')#Change the tick pointing inwards and set the color to blue
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, xlocs=np.arange(100, 123, 4), ylocs=np.arange(0, 25, 2),
linewidth=0.25, linestyle='--', color='k', alpha=0.8)#Add grid lines
gl.top_labels, gl.bottom_labels, gl.right_labels, gl.left_labels = False, False, False, False
plt.savefig('scs_elevation1.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1) # Output the map and set the border margin to be tight
plt.show()
3: Mask other areas and only show the floor plan of a specific area.
# -*- coding: utf-8 -*-
# %%
# Importing related function packages
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as feature
import numpy as np
import matplotlib.ticker as ticker
from cartopy import mpl
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.font_manager import FontProperties
from netCDF4 import Dataset
from palettable.cmocean.diverging import Delta_4
from palettable.colorbrewer.sequential import GnBu_9
from paletable.colorbrewer.sequential import Blues_9
from paletable.scientific.diverging import Roma_20
from pylab import *
def reverse_colourmap(cmap, name='my_cmap_r'):
reverse = []
k = []
for key in cmap._segmentdata:
k.append(key)
channel = cmap._segmentdata[key]
data = []
for t in channel:
data.append((1 - t[0], t[2], t[1]))
reverse.append(sorted(data))
LinearL = dict(zip(k, reverse))
my_cmap_r = mpl.colors.LinearSegmentedColormap(name, LinearL)
return my_cmap_r
cmap = Blues_9.mpl_colormap
cmap_r = reverse_colourmap(cmap)
cmap1 = GnBu_9.mpl_colormap
cmap_r1 = reverse_colourmap(cmap1)
cmap2 = Roma_20.mpl_colormap
cmap_r2 = reverse_colourmap(cmap2)
# read data
a = Dataset('D:\pycharm_work\data\etopo2.nc')
print(a)
lon = a.variables['lon'][:]
lat = a.variables['lat'][:]
ele = a.variables['topo'][:,:]
lon1=lon[1:10800:110]
lat1=lat[1:5400:110]
ele1=ele[1:5400:110,1:10800:110]
print(len(lon1))
print(len(lat))
# Picture 3
#Set map global properties
scale = '50m'
plt.rcParams['font.sans-serif'] = ['Times New Roman'] # Set the overall font to Times New Roman
fig = plt.figure(dpi=300, figsize=(3, 2), facecolor='w', edgecolor='blue')#Set a drawing board and return it to fig
ax = fig.add_axes([0.05, 0.08, 0.92, 0.8], projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([100, 123, 0, 25], crs=ccrs.PlateCarree())# Set the display range
land = feature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
facecolor=feature.COLORS['land'])
ax.add_feature(land, facecolor='0.6')
ax.add_feature(feature.COASTLINE.with_scale('50m'), lw=0.3)#Add coastline: keyword lw sets line width; linestyle sets line style
cs = ax.contourf(lon1, lat1, ele1, levels=np.arange(-9000,0,20),extend='both',cmap=cmap_r1, transform=ccrs.PlateCarree())
# ------colorbar settings
cb = plt.colorbar(cs, ax=ax, extend='both', orientation='vertical',ticks=np.linspace(-9000, 0, 10))
cb.set_label('depth', fontsize=4, color='k')#Set the colorbar label font and size
cb.ax.tick_params(labelsize=4, direction='in') #Set the colorbar tick font size.
cf = ax.contour(lon, lat, ele[:, :], levels=np.linspace(-9000, 0, 5), colors='gray', linestyles='-',
linewidths=0.2,transform=ccrs.PlateCarree())
#Add title
ax.set_title('Etopo', fontsize=4)
# Use Formatter to format tick labels
ax.set_xticks(np.arange(100, 123, 5), crs=ccrs.PlateCarree())#Add latitude and longitude
ax.set_xticklabels(np.arange(100, 123, 5), fontsize=4)
ax.set_yticks(np.arange(0, 25, 5), crs=ccrs.PlateCarree())
ax.set_yticklabels(np.arange(0, 25, 5), fontsize=4)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.tick_params(color='k', direction='in')#Change the tick pointing inwards and set the color to blue
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, xlocs=np.arange(100, 123, 5), ylocs=np.arange(0, 25, 5),
linewidth=0.25, linestyle='--', color='k', alpha=0.8)#Add grid lines
gl.top_labels, gl.bottom_labels, gl.right_labels, gl.left_labels = False, False, False, False
plt.savefig('scs_elevation03.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1) # Output the map and set the border margin to be tight
plt.show()