[PYTHON] Visualisierung von Breiten- / Längenkoordinatendaten (unter der Annahme meteorologischer Daten) unter Verwendung von Cartopy und Matplotlib
Was ist Kartusche?
Cartopy ist ein Python-Paket für die Geodatenverarbeitung zum Zeichnen von Karten und zur Durchführung anderer Geodatenanalysen. Cartopys Zuhause Die Installation ist mit Pip oder Conda einfach.
Hier fassen wir zusammen, wie die Breiten- und Längenkoordinatendaten zusammen mit der Karte gezeichnet werden. Zuerst erkläre ich, wie man eine Karte zeichnet und dann Konturlinien und Vektoren zeichnet.
Basic
Die Module, die häufig zum Zeichnen von Breiten- und Längengraddaten verwendet werden, sind wie folgt. Cartopy kann auch Formdateien zeichnen, daher können Sie sie nach Bedarf importieren.
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.util as cutil
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
Zeichnen Sie eine Karte mit verschiedenen Projektionen
Lassen Sie uns einige typische zeichnen Liste der verwendbaren Projektionen
Global
Zylinderprojektion mit gleichem Abstand und Morweide-Projektion (eine Art Projektion mit gleichem Volumen)
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) #Projektion angeben
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()
Nordpol / Südpolzentrum
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs
fig = plt.figure() #Feigenvorbereitung
#Arktisches Zentrum
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=90.0)
ax = fig.add_subplot(1, 2, 1, projection=proj)
#Zeichenbereich(Längengrad Breitengrad)Bezeichnung von
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())
#Schneiden Sie die Figur in einem Kreis um
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 ")
#Südpolzentrum
proj = ccrs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=-90.0)
ax = fig.add_subplot(1, 2, 2, projection=proj)
#Zeichenbereich(Längengrad Breitengrad)Bezeichnung von
ax.set_extent([-180, 180, -90,-30], ccrs.PlateCarree())
#Schneiden Sie die Figur in einem Kreis um
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 Bereich
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs
fig = plt.figure()
'''
Polare Stereokoordinaten zentriert bei 60 Grad nördlicher Breite und 140 Grad östlicher Länge
'''
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() #Zeigen Sie Land und Meer
ax.coastlines(resolution='50m',) #Erhöhen Sie die Auflösung der Küstenlinie
ax.gridlines()
plt.show()
Zeichnen Sie Längen- / Breitengradkoordinatendaten
Die Daten stammen vom Kyoto University Institute of Survival Area (http://database.rish.kyoto-u.ac.jp/arch/ncep/). Zeichnen Sie die NCEP-Reanalyse aus.
Beispiel-Skriptimportliste
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
Zeichnen Sie eine geopotentiale Höhe von 500 hPa. Da die Daten der Längen- und Breitengradkoordinaten jetzt gezeichnet werden, geben Sie die Projektion an, die Sie in der Achsenprojektion zeichnen möchten, und geben Sie beim Aufrufen von contourf transform = ccrs .. PlateCarree () an. ** Geben Sie transform = ccrs.PlateCarree () ** an, wenn Sie eine Projektion zeichnen. Beachten Sie, dass ** ax.set_extent () ** auch ** ccrs.PlateCarree () ** angibt.
#Lesen Sie 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()
#Daten schneiden
data = hgt[0,5,:,:]
#Zeichnungsteil
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)#Geben Sie den Abstand der Konturlinien an
CN = ax.contour(lon,lat,data,transform=ccrs.PlateCarree(),levels=levels)
ax.clabel(CN,fmt='%.0f') #Beschriften Sie die Konturlinie
#
ax.set_global()
ax.coastlines()
ax.set_title("Contour")
plt.show()
Beim Zeichnen von Breiten- / Längengraddaten werden Konturlinien und Schattierungen bei Unterbrechungen unter periodischen Randbedingungen geschnitten (in diesem Fall Punkte zwischen 357,5 Grad und 360 Grad). Durch Hinzufügen der periodischen Randbedingungen ist es möglich, nahtlos zu zeichnen. Wenn sich der Schnitt am Rand der Figur befindet, ist er nicht erkennbar, sodass Sie ihn möglicherweise nicht hinzufügen müssen.
'''
Die Daten sind die gleichen wie oben
'''
fig = plt.figure(figsize=(20,10))
#Arktisches Zentrum
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())#Zeichenbereich(Längengrad Breitengrad)Bezeichnung von
#Schneiden Sie die Figur in einem Kreis um
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_Geben Sie den Pfad an, um ihn kreisförmig zu machen
ax.coastlines()
plt.colorbar(CF, orientation="horizontal")
ax.set_title( "no add cyclic")
'''
Addition des zyklischen Punktes
'''
ax = fig.add_subplot(1, 2, 2, projection=proj)
ax.set_extent([-180, 180, 30, 90], ccrs.PlateCarree()) #Zeichenbereich(Längengrad Breitengrad)Bezeichnung von
#Schneiden Sie die Figur in einem Kreis um
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)
#Fügen Sie zyklische Punkte hinzu
cyclic_data, cyclic_lon = cutil.add_cyclic_point(data, coord=lon)
#Überprüfen Sie, ob es hinzugefügt wurde
print(lon)
print(cyclic_lon)
#
CF = ax.contourf(cyclic_lon,lat,cyclic_data, transform=ccrs.PlateCarree(),
clip_path=(circle, ax.transAxes) ) # clip_Geben Sie den Pfad an, um ihn kreisförmig zu machen
#
plt.colorbar(CF, orientation="horizontal")
ax.coastlines()
ax.set_title( "add cyclic")
plt.show()
Rechts gibt es keine Brüche in der Figur.
Vektor
200 hPa Wind. Windgeschwindigkeit mit Schatten zeichnen. Hier wird auch die Anpassung der Gitterlinien beschrieben.
#Lesen Sie 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) #Windgeschwindigkeit berechnen
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")
#Extrem, um Fehlermeldungen zu vermeiden(90N,90S)Beschließt, nicht zu zeichnen.
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)
#Vektorlegende anzeigen
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()
#
#Gitterlinieneinstellung
#Verwenden Sie Gitter anstelle von Gitterlinie
ax.set_xticks([0, 60, 120, 180, 240, 300, 359.9999999999], crs=ccrs.PlateCarree()) #0W wird nicht ausgegeben, wenn der Längengrad zum Zeichnen des Gitters auf 360 eingestellt ist
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree()) #Geben Sie den Breitengrad an, um das Raster zu zeichnen
lon_formatter = LongitudeFormatter(zero_direction_label=True) #Längengrad
lat_formatter = LatitudeFormatter(number_format='.1f',degree_symbol='') #Breite. Es ist auch möglich, das Format anzugeben
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.grid()
plt.show()
Rationalisierungsfunktion
Stromlinie von 850 hPa. Stream-Plot verwenden.
#Lesen Sie 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