无法使用 Cartopy 将降水数据从 Stereographic 重新投影到 PlateCarree()
Can't re-project precipitation data from Stereographic to PlateCarree() using Cartopy
我正在尝试绘制来自国家气象局的降水数据。但是，数据默认设置为立体投影。我想在 PlateCarree 投影中绘图，但我遇到了一些困难。当我尝试在 Cartopy 中使用 PlateCarree 投影时，它会绘制地图但不会覆盖降水数据。我假设这意味着我没有正确地将数据从立体投影到 PlateCarree。为了正确地重新投影数据，我需要做些什么吗？
这是与立体投影一起使用的代码：
'''
=====================
NWS Precipitation Map
=====================
Plot a 1-day precipitation map using a netCDF file from the National Weather Service.
This opens the data directly in memory using the support in the netCDF library to open
from an existing memory buffer. In addition to CartoPy and Matplotlib, this uses
a custom colortable as well as MetPy's unit support.
"""
###############################
# Imports
from datetime import datetime, timedelta
from urllib.request import urlopen
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
from metpy.units import masked_array, units
from netCDF4 import Dataset
import pandas as pd
###############################
# Download the data from the National Weather Service.
dt = datetime.utcnow() - timedelta(days=1) # This should always be available
url = 'http://water.weather.gov/precip/downloads/{dt:%Y/%m/%d}/nws_precip_1day_'\
'{dt:%Y%m%d}_conus.nc'.format(dt=dt)
data = urlopen(url).read()
nc = Dataset('data', memory=data)
###############################
# Pull the needed information out of the netCDF file
prcpvar = nc.variables['observation']
data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('in')
#data = data * 0.0393
x = nc.variables['x'][:]
y = nc.variables['y'][:]
proj_var = nc.variables[prcpvar.grid_mapping]
#%%
###############################
# Set up the projection information within CartoPy
globe = ccrs.Globe(semimajor_axis=proj_var.earth_radius)
proj = ccrs.Stereographic(central_latitude=90.0,
central_longitude=proj_var.straight_vertical_longitude_from_pole,
true_scale_latitude=proj_var.standard_parallel, globe=globe)
###############################
# Create the figure and plot the data
# create figure and axes instances
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(111, projection=proj)
#ax.set_extent([-75,-85,35,39])
#draw coastlines, state and country boundaries, edge of map.
ax.coastlines(resolution='10m')
ax.add_feature(cfeature.BORDERS.with_scale('10m'), linewidth=1.5)
ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=2.0)
ax.add_feature(USCOUNTIES.with_scale('500k'), edgecolor='black')
# draw filled contours.
clevs = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
# In future MetPy
# norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)
cmap_data = [
"#04e9e7", # 0.01 - 0.10 inches
"#019ff4", # 0.10 - 0.25 inches
"#0300f4", # 0.25 - 0.50 inches
"#02fd02", # 0.50 - 0.75 inches
"#01c501", # 0.75 - 1.00 inches
"#008e00", # 1.00 - 1.50 inches
"#fdf802", # 1.50 - 2.00 inches
"#e5bc00", # 2.00 - 2.50 inches
"#fd9500", # 2.50 - 3.00 inches
"#fd0000", # 3.00 - 4.00 inches
"#d40000", # 4.00 - 5.00 inches
"#bc0000", # 5.00 - 6.00 inches
"#f800fd", # 6.00 - 8.00 inches
"#9854c6", # 8.00 - 10.00 inches
"#fdfdfd" # 10.00+
]
cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)
cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm)
# add colorbar.
cbar = plt.colorbar(cs, orientation='vertical')
cbar.set_label(data.units)
time = nc.creation_time[4:6]+'/'+nc.creation_time[6:8]+'/'+nc.creation_time[0:4]+' '\
+nc.creation_time[9:11] +':'+ nc.creation_time[11:13] + " UTC"
print(time)
ax.set_title('24 hr Precipitation (in)' + '\n for period ending ' + time, fontsize = 16, fontweight = 'bold' )
'''
但是，当我将投影线更改为 PlateCarree 时，我 运行 出现了上述问题。有人对如何重新投影此数据有任何建议吗？
谢谢
只需将转换参数添加到 ax.contourf 方法。
# Imports
from datetime import datetime, timedelta
from urllib.request import urlopen
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
from metpy.units import masked_array, units
from netCDF4 import Dataset
import pandas as pd
###############################
# Download the data from the National Weather Service.
dt = datetime.utcnow() - timedelta(days=1) # This should always be available
url = 'http://water.weather.gov/precip/downloads/{dt:%Y/%m/%d}/nws_precip_1day_'\
'{dt:%Y%m%d}_conus.nc'.format(dt=dt)
data = urlopen(url).read()
nc = Dataset('data', memory=data)
###############################
# Pull the needed information out of the netCDF file
prcpvar = nc.variables['observation']
data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('in')
#data = data * 0.0393
x = nc.variables['x'][:]
y = nc.variables['y'][:]
proj_var = nc.variables[prcpvar.grid_mapping]
#%%
###############################
# Set up the projection information within CartoPy
globe = ccrs.Globe(semimajor_axis=proj_var.earth_radius)
proj = ccrs.Stereographic(central_latitude=90.0,
central_longitude=proj_var.straight_vertical_longitude_from_pole,
true_scale_latitude=proj_var.standard_parallel, globe=globe)
###############################
# Create the figure and plot the data
# create figure and axes instances
fig = plt.figure(figsize=(15, 15))
pc_proj = ccrs.PlateCarree()
#ax = fig.add_subplot(111, projection=proj)
ax = fig.add_subplot(111, projection=pc_proj)
ax.set_extent([-128,-65,25,52])
#draw coastlines, state and country boundaries, edge of map.
ax.coastlines(resolution='10m')
#ax.add_feature(cfeature.BORDERS.with_scale('10m'), linewidth=1.5)
#ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=2.0)
#ax.add_feature(USCOUNTIES.with_scale('500k'), edgecolor='black')
# draw filled contours.
clevs = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
# In future MetPy
# norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)
cmap_data = [
"#04e9e7", # 0.01 - 0.10 inches
"#019ff4", # 0.10 - 0.25 inches
"#0300f4", # 0.25 - 0.50 inches
"#02fd02", # 0.50 - 0.75 inches
"#01c501", # 0.75 - 1.00 inches
"#008e00", # 1.00 - 1.50 inches
"#fdf802", # 1.50 - 2.00 inches
"#e5bc00", # 2.00 - 2.50 inches
"#fd9500", # 2.50 - 3.00 inches
"#fd0000", # 3.00 - 4.00 inches
"#d40000", # 4.00 - 5.00 inches
"#bc0000", # 5.00 - 6.00 inches
"#f800fd", # 6.00 - 8.00 inches
"#9854c6", # 8.00 - 10.00 inches
"#fdfdfd" # 10.00+
]
cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)
#cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm)
# add transform args
cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm,transform=proj)
# add colorbar.
cbar = plt.colorbar(cs, orientation='vertical')
cbar.set_label(data.units)
time = nc.creation_time[4:6]+'/'+nc.creation_time[6:8]+'/'+nc.creation_time[0:4]+' '\
+nc.creation_time[9:11] +':'+ nc.creation_time[11:13] + " UTC"
print(time)
ax.set_title('24 hr Precipitation (in)' + '\n for period ending ' + time, fontsize = 16, fontweight = 'bold' )
plt.savefig("prec_usa_pc.png")
#plt.show()
下面是输出图
我正在尝试绘制来自国家气象局的降水数据。但是，数据默认设置为立体投影。我想在 PlateCarree 投影中绘图，但我遇到了一些困难。当我尝试在 Cartopy 中使用 PlateCarree 投影时，它会绘制地图但不会覆盖降水数据。我假设这意味着我没有正确地将数据从立体投影到 PlateCarree。为了正确地重新投影数据，我需要做些什么吗？
这是与立体投影一起使用的代码：
'''
=====================
NWS Precipitation Map
=====================
Plot a 1-day precipitation map using a netCDF file from the National Weather Service.
This opens the data directly in memory using the support in the netCDF library to open
from an existing memory buffer. In addition to CartoPy and Matplotlib, this uses
a custom colortable as well as MetPy's unit support.
"""
###############################
# Imports
from datetime import datetime, timedelta
from urllib.request import urlopen
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
from metpy.units import masked_array, units
from netCDF4 import Dataset
import pandas as pd
###############################
# Download the data from the National Weather Service.
dt = datetime.utcnow() - timedelta(days=1) # This should always be available
url = 'http://water.weather.gov/precip/downloads/{dt:%Y/%m/%d}/nws_precip_1day_'\
'{dt:%Y%m%d}_conus.nc'.format(dt=dt)
data = urlopen(url).read()
nc = Dataset('data', memory=data)
###############################
# Pull the needed information out of the netCDF file
prcpvar = nc.variables['observation']
data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('in')
#data = data * 0.0393
x = nc.variables['x'][:]
y = nc.variables['y'][:]
proj_var = nc.variables[prcpvar.grid_mapping]
#%%
###############################
# Set up the projection information within CartoPy
globe = ccrs.Globe(semimajor_axis=proj_var.earth_radius)
proj = ccrs.Stereographic(central_latitude=90.0,
central_longitude=proj_var.straight_vertical_longitude_from_pole,
true_scale_latitude=proj_var.standard_parallel, globe=globe)
###############################
# Create the figure and plot the data
# create figure and axes instances
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(111, projection=proj)
#ax.set_extent([-75,-85,35,39])
#draw coastlines, state and country boundaries, edge of map.
ax.coastlines(resolution='10m')
ax.add_feature(cfeature.BORDERS.with_scale('10m'), linewidth=1.5)
ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=2.0)
ax.add_feature(USCOUNTIES.with_scale('500k'), edgecolor='black')
# draw filled contours.
clevs = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
# In future MetPy
# norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)
cmap_data = [
"#04e9e7", # 0.01 - 0.10 inches
"#019ff4", # 0.10 - 0.25 inches
"#0300f4", # 0.25 - 0.50 inches
"#02fd02", # 0.50 - 0.75 inches
"#01c501", # 0.75 - 1.00 inches
"#008e00", # 1.00 - 1.50 inches
"#fdf802", # 1.50 - 2.00 inches
"#e5bc00", # 2.00 - 2.50 inches
"#fd9500", # 2.50 - 3.00 inches
"#fd0000", # 3.00 - 4.00 inches
"#d40000", # 4.00 - 5.00 inches
"#bc0000", # 5.00 - 6.00 inches
"#f800fd", # 6.00 - 8.00 inches
"#9854c6", # 8.00 - 10.00 inches
"#fdfdfd" # 10.00+
]
cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)
cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm)
# add colorbar.
cbar = plt.colorbar(cs, orientation='vertical')
cbar.set_label(data.units)
time = nc.creation_time[4:6]+'/'+nc.creation_time[6:8]+'/'+nc.creation_time[0:4]+' '\
+nc.creation_time[9:11] +':'+ nc.creation_time[11:13] + " UTC"
print(time)
ax.set_title('24 hr Precipitation (in)' + '\n for period ending ' + time, fontsize = 16, fontweight = 'bold' )
'''
但是，当我将投影线更改为 PlateCarree 时，我 运行 出现了上述问题。有人对如何重新投影此数据有任何建议吗？
谢谢
只需将转换参数添加到 ax.contourf 方法。
# Imports
from datetime import datetime, timedelta
from urllib.request import urlopen
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
from metpy.units import masked_array, units
from netCDF4 import Dataset
import pandas as pd
###############################
# Download the data from the National Weather Service.
dt = datetime.utcnow() - timedelta(days=1) # This should always be available
url = 'http://water.weather.gov/precip/downloads/{dt:%Y/%m/%d}/nws_precip_1day_'\
'{dt:%Y%m%d}_conus.nc'.format(dt=dt)
data = urlopen(url).read()
nc = Dataset('data', memory=data)
###############################
# Pull the needed information out of the netCDF file
prcpvar = nc.variables['observation']
data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('in')
#data = data * 0.0393
x = nc.variables['x'][:]
y = nc.variables['y'][:]
proj_var = nc.variables[prcpvar.grid_mapping]
#%%
###############################
# Set up the projection information within CartoPy
globe = ccrs.Globe(semimajor_axis=proj_var.earth_radius)
proj = ccrs.Stereographic(central_latitude=90.0,
central_longitude=proj_var.straight_vertical_longitude_from_pole,
true_scale_latitude=proj_var.standard_parallel, globe=globe)
###############################
# Create the figure and plot the data
# create figure and axes instances
fig = plt.figure(figsize=(15, 15))
pc_proj = ccrs.PlateCarree()
#ax = fig.add_subplot(111, projection=proj)
ax = fig.add_subplot(111, projection=pc_proj)
ax.set_extent([-128,-65,25,52])
#draw coastlines, state and country boundaries, edge of map.
ax.coastlines(resolution='10m')
#ax.add_feature(cfeature.BORDERS.with_scale('10m'), linewidth=1.5)
#ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=2.0)
#ax.add_feature(USCOUNTIES.with_scale('500k'), edgecolor='black')
# draw filled contours.
clevs = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
# In future MetPy
# norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)
cmap_data = [
"#04e9e7", # 0.01 - 0.10 inches
"#019ff4", # 0.10 - 0.25 inches
"#0300f4", # 0.25 - 0.50 inches
"#02fd02", # 0.50 - 0.75 inches
"#01c501", # 0.75 - 1.00 inches
"#008e00", # 1.00 - 1.50 inches
"#fdf802", # 1.50 - 2.00 inches
"#e5bc00", # 2.00 - 2.50 inches
"#fd9500", # 2.50 - 3.00 inches
"#fd0000", # 3.00 - 4.00 inches
"#d40000", # 4.00 - 5.00 inches
"#bc0000", # 5.00 - 6.00 inches
"#f800fd", # 6.00 - 8.00 inches
"#9854c6", # 8.00 - 10.00 inches
"#fdfdfd" # 10.00+
]
cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)
#cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm)
# add transform args
cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm,transform=proj)
# add colorbar.
cbar = plt.colorbar(cs, orientation='vertical')
cbar.set_label(data.units)
time = nc.creation_time[4:6]+'/'+nc.creation_time[6:8]+'/'+nc.creation_time[0:4]+' '\
+nc.creation_time[9:11] +':'+ nc.creation_time[11:13] + " UTC"
print(time)
ax.set_title('24 hr Precipitation (in)' + '\n for period ending ' + time, fontsize = 16, fontweight = 'bold' )
plt.savefig("prec_usa_pc.png")
#plt.show()
下面是输出图