[PYTHON] Visualization of latitude / longitude coordinate data (assuming meteorological data) using cartopy and matplotlib
What is cartopy
Cartopy is a Python package for geospatial data processing for drawing maps and performing other geospatial data analysis. Cartopy's home Installation is easy with pip or conda.
Here, we will summarize how to draw latitude / longitude coordinate data together with a map. First, I will explain how to draw a map, and then draw contour lines and vectors.
Basic
The modules that are often used to draw latitude / longitude data are as follows. cartopy can also draw shape files, so I think you can import it as needed.
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.util as cutil
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
Draw a map with various projections
Let's plot some typical ones List of usable projections
Global
Equirectangular projection and Mollweide projection (a type of equal volume projection)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig = plt.figure(figsize=(10,20))
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 2, 1, projection=proj) #Specify projection
ax.set_global()
ax.coastlines()
ax.gridlines()
ax.set_title("PlateCarree")
# mollweide
proj = ccrs.Mollweide()
ax = fig.add_subplot(1, 2, 2, projection=proj)
ax.set_global()
ax.coastlines()
ax.gridlines()
ax.set_title("mollweide")
plt.show()
Arctic / Antarctic center
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs
fig = plt.figure() #fig preparation
#Arctic center
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=90.0)
ax = fig.add_subplot(1, 2, 1, projection=proj)
#Drawing range(longitude latitude)Designation of
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())
#Cut the circumference of the figure into a circle
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
#
ax.coastlines()
ax.gridlines()
ax.set_title( " NP ")
#Antarctica center
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=-90.0)
ax = fig.add_subplot(1, 2, 2, projection=proj)
#Drawing range(longitude latitude)Designation of
ax.set_extent([-180, 180, -90,-30], ccrs.PlateCarree())
#Cut the circumference of the figure into a circle
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
#
ax.coastlines()
ax.gridlines()
ax.set_title( " SP ")
plt.show()
Japan area
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs
fig = plt.figure()
'''
Polar stereo coordinates centered at latitude 60 degrees north and longitude 140 degrees east
'''
proj = ccrs.Stereographic(central_latitude=60, central_longitude=140)
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.set_extent([120, 160, 20, 50], ccrs.PlateCarree())
#ax.set_global()
ax.stock_img() #Show land and sea
ax.coastlines(resolution='50m',) #Increase coastline resolution
ax.gridlines()
plt.show()
Draw latitude / longitude coordinate data
Data is from the Research Institute for Sustainable Humanosphere, Kyoto University (http://database.rish.kyoto-u.ac.jp/arch/ncep/) Draw the NCEP reanalysis obtained from.
Sample script import list
import netCDF4
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import matplotlib.cm as cm
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.util as cutil
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
Isoline
Draws a geopotential height of 500hPa. Since the data of latitude and longitude coordinates is drawn now, specify the projection you want to draw in the axes projection, and specify transform = ccrs .. PlateCarree () when calling contourf. ** Specify transform = ccrs.PlateCarree () ** when drawing any projection. Note that ** ax.set_extent () ** also specifies ** ccrs.PlateCarree () **.
#read netcdf
nc = netCDF4.Dataset("hgt.mon.ltm.nc","r")
hgt = nc.variables["hgt"][:]
level = nc.variables["level"][:]
lat = nc.variables["lat"][:]
lon = nc.variables["lon"][:]
nc.close()
#Data cut
data = hgt[0,5,:,:]
#Drawing part
fig = plt.figure(figsize=(10,10))
proj = ccrs.PlateCarree(central_longitude= 180)
ax = fig.add_subplot(1, 1, 1, projection=proj)
#
levels=np.arange(5100,5800,60)#Specify the spacing of contour lines
CN = ax.contour(lon,lat,data,transform=ccrs.PlateCarree(),levels=levels)
ax.clabel(CN,fmt='%.0f') #Label the contour line
#
ax.set_global()
ax.coastlines()
ax.set_title("Contour")
plt.show()
When drawing latitude / longitude data, contour lines and shades are cut at the breaks in the periodic boundary conditions (in this case, points between 357.5 degrees and 360 degrees). By adding the periodic boundary conditions, it is possible to draw seamlessly. If the cut is at the edge of the figure, it will not be noticeable, so you may not need to add it.
'''
Data is the same as above
'''
fig = plt.figure(figsize=(20,10))
#Arctic center
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=90.0)
ax = fig.add_subplot(1, 2, 1, projection=proj)
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())#Drawing range(longitude latitude)Designation of
#Cut the circumference of the figure into a circle
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
#
CF = ax.contourf(lon,lat,hgt[0,5,:,:], transform=ccrs.PlateCarree(),
clip_path=(circle, ax.transAxes) ) # clip_Specify path to make it circular
ax.coastlines()
plt.colorbar(CF, orientation="horizontal")
ax.set_title( "no add cyclic")
'''
Addition of cyclic point
'''
ax = fig.add_subplot(1, 2, 2, projection=proj)
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree()) #Drawing range(longitude latitude)Designation of
#Cut the circumference of the figure into a circle
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
#Add cyclic points
cyclic_data, cyclic_lon = cutil.add_cyclic_point(data, coord=lon)
#Check if it was added
print(lon)
print(cyclic_lon)
#
CF = ax.contourf(cyclic_lon,lat,cyclic_data, transform=ccrs.PlateCarree(),
clip_path=(circle, ax.transAxes) ) # clip_Specify path to make it circular
#
plt.colorbar(CF, orientation="horizontal")
ax.coastlines()
ax.set_title( "add cyclic")
plt.show()
On the right, there are no breaks in the figure.
vector
200 hPa wind. Draw the wind speed with shade. The adjustment of grid lines is also described here.
#read netcdf
nc = netCDF4.Dataset("uwnd.mon.ltm.nc","r")
u = nc.variables["uwnd"][:][0,9,:,:]
level = nc.variables["level"][:]
lat = nc.variables["lat"][:]
lon = nc.variables["lon"][:]
nc.close()
#
nc = netCDF4.Dataset("vwnd.mon.ltm.nc","r")
v = nc.variables["vwnd"][:][0,9,:,:]
level = nc.variables["level"][:]
lat = nc.variables["lat"][:]
lon = nc.variables["lon"][:]
nc.close()
#
fig = plt.figure(figsize=(10,10))
proj = ccrs.PlateCarree(central_longitude= 180)
ax = fig.add_subplot(1, 1, 1, projection=proj)
#
sp = np.sqrt(u**2+v**2) #Calculate wind speed
sp, cyclic_lon = cutil.add_cyclic_point(sp, coord=lon)
#
levels=np.arange(0,61,5)
cf = ax.contourf(cyclic_lon, lat, sp, transform=ccrs.PlateCarree(), levels=levels, cmap=cm.jet, extend = "both")
#Extreme to avoid error messages(90N,90S)Decides not to draw.
Q = ax.quiver(lon,lat[1:-1],u[1:-1,:],v[1:-1,:],transform=ccrs.PlateCarree() , regrid_shape=20, units='xy', angles='xy', scale_units='xy', scale=1)
#Display vector legend
qk = ax.quiverkey(Q, 0.8, 1.05, 20, r'$20 m/s$', labelpos='E',
coordinates='axes',transform=ccrs.PlateCarree() )
plt.colorbar(cf, orientation="horizontal" )
#
ax.coastlines()
ax.set_title("Vector and Wind speed(shade)")
ax.set_global()
#
#Grid line adjustment
#Use grid instead of gridline
ax.set_xticks([0, 60, 120, 180, 240, 300, 359.9999999999], crs=ccrs.PlateCarree()) #0W does not come out when the longitude to draw grid is set to 360
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree()) #Specify the latitude to draw the grid
lon_formatter = LongitudeFormatter(zero_direction_label=True) #longitude
lat_formatter = LatitudeFormatter(number_format='.1f',degree_symbol='') #latitude. It is also possible to specify format
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.grid()
plt.show()
Streamline function
Streamline of 850 hPa. Use stream plot.
#read netcdf
nc = netCDF4.Dataset("uwnd.mon.ltm.nc","r")
u = nc.variables["uwnd"][:][7,2,:,:]
level = nc.variables["level"][:]
lat = nc.variables["lat"][:]
lon = nc.variables["lon"][:]
nc.close()
#
nc = netCDF4.Dataset("vwnd.mon.ltm.nc","r")
v = nc.variables["vwnd"][:][7,2,:,:]
level = nc.variables["level"][:]
lat = nc.variables["lat"][:]
lon = nc.variables["lon"][:]
nc.close()
'''
plot
'''
fig = plt.figure(figsize=(10,10))
proj = ccrs.PlateCarree(central_longitude= 180)
ax = fig.add_subplot(1, 1, 1, projection=proj)
stream = ax.streamplot(lon,lat,u,v,transform=ccrs.PlateCarree())
ax.clabel(CN,fmt='%.0f')
ax.set_global()
ax.coastlines()
ax.set_title("Stream line")
plt.show()
Recommended Posts