[PYTHON] Visualisation des données de coordonnées de latitude / longitude (en supposant des données météorologiques) à l'aide de cartopy et matplotlib

Qu'est-ce que la cartopie

Cartopy est un package Python pour le traitement des données géospatiales pour dessiner des cartes et effectuer d'autres analyses de données géospatiales. Accueil de Cartopy L'installation est facile avec pip ou conda.

Ici, nous allons résumer comment dessiner les données de coordonnées de latitude et de longitude avec la carte. Tout d'abord, je vais expliquer comment dessiner une carte, puis dessiner des courbes de niveau et des vecteurs.

De base

Les modules qui sont souvent utilisés pour dessiner des données de latitude et de longitude sont les suivants. cartopy peut également dessiner des fichiers de forme, donc je pense que vous pouvez les importer selon vos besoins.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.util as cutil
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

Dessinez une carte avec diverses projections

Tracons quelques exemples typiques Liste des projections utilisables

Global

Projection cylindrique à égale distance et projection Morweide (un type de projection à volume égal)

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) #Spécifiez la 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()

Pôle Nord / Pôle Sud Centre

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs

fig = plt.figure() #préparation de figues

#Centre arctique
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=90.0)
ax = fig.add_subplot(1, 2, 1, projection=proj) 
#Plage de dessin(longitude latitude)Désignation de
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())
#Couper autour de la figure en cercle
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 ")

#Centre du pôle Sud
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=-90.0)
ax = fig.add_subplot(1, 2, 2, projection=proj) 
#Plage de dessin(longitude latitude)Désignation de
ax.set_extent([-180, 180, -90,-30], ccrs.PlateCarree())
#Couper autour de la figure en cercle
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()

Région du Japon

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs

fig = plt.figure()
'''
Coordonnées stéréo polaires centrées à 60 degrés de latitude nord et 140 degrés de longitude est
'''
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() #Afficher la terre et la mer
ax.coastlines(resolution='50m',) #Augmenter la résolution du littoral
ax.gridlines()
plt.show()

Dessiner des données de coordonnées de latitude / longitude

Les données proviennent de l'Institut de survie de l'Université de Kyoto (http://database.rish.kyoto-u.ac.jp/arch/ncep/) Dessinez la réanalyse NCEP obtenue à partir de.

Exemple de liste d'importation de script

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

Dessinez une altitude géopotentielle de 500 hPa. Puisque les données de coordonnées de latitude et de longitude sont maintenant dessinées, spécifiez la projection que vous voulez dessiner dans la projection des axes et spécifiez transform = ccrs..PlateCarree () lors de l'appel de contourf. ** Spécifiez transform = ccrs.PlateCarree () ** lorsque vous dessinez une projection. Notez que ** ax.set_extent () ** spécifie également ** ccrs.PlateCarree () **.

#Lire 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()
#Couper les données
data = hgt[0,5,:,:]
#Partie de dessin
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)#Spécifiez l'espacement des courbes de niveau
CN = ax.contour(lon,lat,data,transform=ccrs.PlateCarree(),levels=levels)
ax.clabel(CN,fmt='%.0f') #Étiqueter la ligne de contour
#
ax.set_global()
ax.coastlines()
ax.set_title("Contour")
plt.show()

Lorsque vous dessinez des données de latitude / longitude, les lignes de contour et les ombres sont coupées aux ruptures dans des conditions aux limites périodiques (dans ce cas, des points entre 357,5 degrés et 360 degrés). En ajoutant les conditions aux limites périodiques, il est possible de dessiner de manière transparente. Si la coupe est au bord de la figure, elle ne sera pas perceptible, vous n'aurez donc peut-être pas besoin de l'ajouter.

'''
Les données sont les mêmes que ci-dessus
'''
fig = plt.figure(figsize=(20,10))
#Centre arctique
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())#Plage de dessin(longitude latitude)Désignation de
#Couper autour de la figure en cercle
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_Spécifiez le chemin pour le rendre circulaire
ax.coastlines()
plt.colorbar(CF, orientation="horizontal")
ax.set_title( "no add cyclic")
'''
Ajout de point cyclique
'''
ax = fig.add_subplot(1, 2, 2, projection=proj) 
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())　#Plage de dessin(longitude latitude)Désignation de
#Couper autour de la figure en cercle
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)
#Ajouter des points cycliques
cyclic_data, cyclic_lon = cutil.add_cyclic_point(data, coord=lon)
#Vérifiez s'il a été ajouté
print(lon)
print(cyclic_lon)
# 
CF = ax.contourf(cyclic_lon,lat,cyclic_data, transform=ccrs.PlateCarree(),
                clip_path=(circle, ax.transAxes) ) # clip_Spécifiez le chemin pour le rendre circulaire
#
plt.colorbar(CF, orientation="horizontal")
ax.coastlines()
ax.set_title( "add cyclic")
plt.show()

Sur la droite, il n'y a pas de coupures dans la figure.

vecteur

Vent de 200 hPa. Dessinez la vitesse du vent avec de l'ombre. L'ajustement des lignes de quadrillage est également décrit ici.

#Lire 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) #Calculer la vitesse du vent
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")
#Extrême pour éviter les messages d'erreur(90N,90S)Décide de ne pas dessiner.
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)
#Afficher la légende du vecteur
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()
#
#Réglage de la ligne de grille
#Utiliser la grille au lieu du quadrillage
ax.set_xticks([0, 60, 120, 180, 240, 300, 359.9999999999], crs=ccrs.PlateCarree()) #0W ne sort pas lorsque la longitude pour dessiner la grille est définie sur 360
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree()) #Spécifiez la latitude pour dessiner la grille
lon_formatter = LongitudeFormatter(zero_direction_label=True) #longitude
lat_formatter = LatitudeFormatter(number_format='.1f',degree_symbol='') #latitude. Il est également possible de spécifier le format
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.grid()
plt.show()

Fonction de rationalisation

Rationalisation de 850 hPa. Utilisez le tracé de flux.

#Lire 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()