英文:
How do I plot netCDF sea ice extent dataset in Python?
问题
我计划使用Python绘制海洋冰盖在南极附近的分布图，使用NOAA-NSIDC被动微波海冰浓度第4版本数据集。我正在使用以下代码：
import xarray as xr
import numpy as np
import cmocean
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
# 读取netCDF文件
ds = xr.open_dataset('seaice_conc_monthly_sh_202212_f17_v04r00.nc')
# 提取第一个时间步的海冰范围数据
sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values
# 提取x和y网格坐标
xgrid = ds.xgrid.values
ygrid = ds.ygrid.values
# 创建图形
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=ccrs.PlateCarree(), cmap='cmo.ice')
ax.gridlines()
ax.coastlines()
gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
cbar = fig.colorbar(im, ax=ax, shrink=0.4)
cbar.ax.set_ylabel('Sea Ice Concentration (%)', fontsize=16)
plt.show()
然而，我遇到了以下错误：
TypeError: 'GeometryCollection' object is not subscriptable
请问如何在适合南极的投影下绘制此图？以下是有关数据集的一些信息：
数据集中包含以下变量：
- cdr_seaice_conc_monthly
- nsidc_bt_seaice_conc_monthly
- nsidc_nt_seaice_conc_monthly
- projection
- qa_of_cdr_seaice_conc_monthly
- stdev_of_cdr_seaice_conc_monthly
我想要获得类似于此图的绘图：
英文:
I am planning to plot sea ice extent around Antarctica using NOAA-NSIDC Passive Microwave Sea Ice Concentration, Version 4 dataset in Python. I am using the following code:
>>> import xarray as xr
>>> import numpy as np
>>> import cmocean
>>> import cartopy.crs as ccrs
>>> import matplotlib.pyplot as plt
>>> import cartopy.feature as cfeature
>>> from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
>>> # Read in the netCDF file
>>>
>>> ds = xr.open_dataset('seaice_conc_monthly_sh_202212_f17_v04r00.nc')
>>> # Extract the sea ice extent data for the first time step
>>>
>>> sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values
>>> # Extract the x and y grid coordinates
>>>
>>> xgrid = ds.xgrid.values
>>> ygrid = ds.ygrid.values
>>> # Create the plot
>>>
>>> fig = plt.figure(figsize=(10, 10))
>>> ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
>>> ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
>>> im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=ccrs.PlateCarree(), cmap='cmo.ice')
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 318, in wrapper
return func(self, *args, **kwargs)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 1787, in pcolormesh
result = self._wrap_quadmesh(result, **kwargs)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 1970, in _wrap_quadmesh
pcolor_col = self.pcolor(coords[..., 0], coords[..., 1],
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 318, in wrapper
return func(self, *args, **kwargs)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 2015, in pcolor
limits = result.get_datalim(self.transData)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 266, in get_datalim
paths = [transform.transform_path_non_affine(p) for p in paths]
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 266, in <listcomp>
paths = [transform.transform_path_non_affine(p) for p in paths]
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\transforms.py", line 2431, in transform_path_non_affine
return self._a.transform_path_non_affine(path)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 190, in transform_path_non_affine
proj_geom = self.target_projection.project_geometry(
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 808, in project_geometry
return getattr(self, method_name)(geometry, src_crs)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 963, in _project_polygon
return self._rings_to_multi_polygon(rings, is_ccw)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 1224, in _rings_to_multi_polygon
multi_poly = sgeom.MultiPolygon(polygon_bits)
File "C:\Users\ashif\miniconda3\lib\site-packages\shapely\geometry\multipolygon.py", line 79, in __new__
shell = ob[0]
TypeError: 'GeometryCollection' object is not subscriptable
>>> ax.gridlines()
<cartopy.mpl.gridliner.Gridliner object at 0x000001EE3C35B8E0>
>>> ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist object at 0x000001EE3C31FE80>
>>> gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
>>> gl.top_labels = False
>>> gl.right_labels = False
>>> gl.xformatter = LONGITUDE_FORMATTER
>>> gl.yformatter = LATITUDE_FORMATTER
>>> cbar = fig.colorbar(im, ax=ax, shrink=0.4)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
NameError: name 'im' is not defined. Did you mean: 'id'?
>>> cbar.ax.set_ylabel('Sea Ice Concentration (%)', fontsize=16)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
NameError: name 'cbar' is not defined
>>> plt.show()
Traceback (most recent call last):
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\backends\backend_qt.py", line 468, in _draw_idle
self.draw()
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\backends\backend_agg.py", line 400, in draw
self.figure.draw(self.renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 95, in draw_wrapper
result = draw(artist, renderer, *args, **kwargs)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
return draw(artist, renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\figure.py", line 3140, in draw
mimage._draw_list_compositing_images(
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\image.py", line 131, in _draw_list_compositing_images
a.draw(renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
return draw(artist, renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 538, in draw
return super().draw(renderer=renderer, **kwargs)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
return draw(artist, renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\axes\_base.py", line 3064, in draw
mimage._draw_list_compositing_images(
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\image.py", line 131, in _draw_list_compositing_images
a.draw(renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
return draw(artist, renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 972, in draw
super().draw(renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
return draw(artist, renderer)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 351, in draw
transform, offset_trf, offsets, paths = self._prepare_points()
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 328, in _prepare_points
paths = [transform.transform_path_non_affine(path)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 328, in <listcomp>
paths = [transform.transform_path_non_affine(path)
File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\transforms.py", line 2431, in transform_path_non_affine
return self._a.transform_path_non_affine(path)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 190, in transform_path_non_affine
proj_geom = self.target_projection.project_geometry(
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 808, in project_geometry
return getattr(self, method_name)(geometry, src_crs)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 963, in _project_polygon
return self._rings_to_multi_polygon(rings, is_ccw)
File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 1224, in _rings_to_multi_polygon
multi_poly = sgeom.MultiPolygon(polygon_bits)
File "C:\Users\ashif\miniconda3\lib\site-packages\shapely\geometry\multipolygon.py", line 79, in __new__
shell = ob[0]
TypeError: 'GeometryCollection' object is not subscriptable
However I am facing the following error:
TypeError: 'GeometryCollection' object is not subscriptable
Could anyone tell me how do I plot this in a projection that is suitable for the south poles? Below are some of the information about the dataset:
Attributes:
long_name: NOAA/NSIDC Climate Data Record of Passive Microwave...
standard_name: sea_ice_area_fraction
units: 1
flag_values: [251 252 253 254 255]
flag_meanings: pole_hole lakes coastal land_mask missing_data
datum: +ellps=urn:ogc:def:crs:EPSG::4326
grid_mapping: projection
reference: https://nsidc.org/data/g02202/versions/4/
ancillary_variables: stdev_of_cdr_seaice_conc_monthly qa_of_cdr_seaice_c...
valid_range: [ 0 100], 'nsidc_bt_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
[104912 values with dtype=float32]
Attributes:
long_name: Passive Microwave Monthly Southern Hemisphere Sea Ice Con...
standard_name: sea_ice_area_fraction
units: 1
flag_values: [251 252 253 254 255]
flag_meanings: pole_hole unused coastal land_mask missing_data
datum: +ellps=urn:ogc:def:crs:EPSG::4326
grid_mapping: projection
valid_range: [ 0 100], 'nsidc_nt_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
[104912 values with dtype=float32]
Attributes:
long_name: Passive Microwave Monthly Southern Hemisphere Sea Ice Con...
standard_name: sea_ice_area_fraction
units: 1
flag_values: [251 252 253 254 255]
flag_meanings: pole_hole unused coastal land_mask missing_data
datum: +ellps=urn:ogc:def:crs:EPSG::4326
grid_mapping: projection
valid_range: [ 0 100], 'projection': <xarray.Variable ()>
[1 values with dtype=|S1]
Attributes: (12/23)
grid_boundary_top_projected_y: 4350000.0
grid_boundary_bottom_projected_y: -3950000.0
grid_boundary_right_projected_x: 3950000.0
grid_boundary_left_projected_x: -3950000.0
parent_grid_cell_row_subset_start: 0.0
parent_grid_cell_row_subset_end: 332.0
... ...
scaling_factor: 1.0
false_easting: 0.0
false_northing: 0.0
semimajor_radius: 6378273.0
semiminor_radius: 6356889.449
units: meters, 'qa_of_cdr_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
[104912 values with dtype=float32]
Attributes:
long_name: Passive Microwave Monthly Southern Hemisphere Sea Ice Con...
standard_name: sea_ice_area_fraction status_flag
flag_meanings: Average_concentration_exceeds_0.15 Average_concentration_...
datum: +ellps=urn:ogc:def:crs:EPSG::4326
grid_mapping: projection
flag_masks: [ 1 2 4 8 32 64 128]
valid_range: [ 1 255], 'stdev_of_cdr_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
[104912 values with dtype=float32]
Attributes:
long_name: Passive Microwave Monthly Southern Hemisphere Sea Ice Conc...
datum: +ellps=urn:ogc:def:crs:EPSG::4326
grid_mapping: projection
valid_range: [0. 1.], 'time': <xarray.Variable (tdim: 1)>
[1 values with dtype=datetime64[ns]]
Attributes:
standard_name: time
long_name: ANSI date
axis: T, 'xgrid': <xarray.Variable (x: 316)>
array([-3937500., -3912500., -3887500., ..., 3887500., 3912500., 3937500.],
dtype=float32)
Attributes:
valid_range: [-3950000. 3950000.]
units: meters
long_name: projection_grid_x_centers
standard_name: projection_x_coordinate
axis: X, 'ygrid': <xarray.Variable (y: 332)>
array([ 4337500., 4312500., 4287500., ..., -3887500., -3912500., -3937500.],
dtype=float32)
Attributes:
valid_range: [-3950000. 4350000.]
units: meters
long_name: projection_grid_y_centers
standard_name: projection_y_coordinate
axis: Y}))
The following variables are in the dataset:
cdr_seaice_conc_monthly (tdim, y, x) float32 ...
nsidc_bt_seaice_conc_monthly (tdim, y, x) float32 ...
nsidc_nt_seaice_conc_monthly (tdim, y, x) float32 ...
projection |S1 ...
qa_of_cdr_seaice_conc_monthly (tdim, y, x) float32 ...
stdev_of_cdr_seaice_conc_monthly (tdim, y, x) float32 ...
I would like to get a plot similar to this:
答案1
得分: 0
主要问题在于您指定了数据的投影为ccrs.PlateCarree()，这显然是不正确的。
也许有比下面示例中解析更好的方法，但它捕捉到了要点：
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.patheffects as PathEffects
with xr.open_dataset("seaice_conc_monthly_sh_202212_f17_v04r00.nc") as ds:
sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values * 100
data_proj = ccrs.epsg(ds.projection.attrs["srid"].split("::")[-1])
flags = ds.cdr_seaice_conc_monthly.attrs["flag_values"]
invalid = np.isin(sea_ice_extent, flags)
sea_ice_extent[invalid] = np.nan
xgrid = ds.xgrid.values
ygrid = ds.ygrid.values
cmap = mpl.colormaps["BuPu"].copy()
cmap.set_bad("grey", alpha=.2)
fig, ax = plt.subplots(
figsize=(8,8), facecolor="w",
subplot_kw=dict(projection=ccrs.SouthPolarStereo()),
)
im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=data_proj, cmap=cmap)
ax.set_extent([-180, 180, -90, -45], ccrs.PlateCarree())
ax.coastlines()
gl = ax.gridlines(
draw_labels=True,
lw=0.4,
ylabel_style=dict(path_effects=[PathEffects.withStroke(linewidth=2, foreground="w")]),
xlocs=range(-180,180,45), ylocs=range(-90,0,10),
y_inline=True,
)
cbar = fig.colorbar(im, ax=ax, shrink=0.4, label='Sea Ice Concentration (%)')
ax.set_boundary(mpl.path.Path.circle((0, 0), 3950000*1.2), transform=data_proj)
英文:
The main issue is that you're specifying the projection of the data to be ccrs.PlateCarree(), which clearly is incorrect.
There might be a better way than to parse it as in the example below, but it captures the gist:
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.patheffects as PathEffects
with xr.open_dataset("seaice_conc_monthly_sh_202212_f17_v04r00.nc") as ds:
sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values * 100
data_proj = ccrs.epsg(ds.projection.attrs["srid"].split("::")[-1])
flags = ds.cdr_seaice_conc_monthly.attrs["flag_values"]
invalid = np.isin(sea_ice_extent, flags)
sea_ice_extent[invalid] = np.nan
xgrid = ds.xgrid.values
ygrid = ds.ygrid.values
cmap = mpl.colormaps["BuPu"].copy()
cmap.set_bad("grey", alpha=.2)
fig, ax = plt.subplots(
figsize=(8,8), facecolor="w",
subplot_kw=dict(projection=ccrs.SouthPolarStereo()),
)
im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=data_proj, cmap=cmap)
ax.set_extent([-180, 180, -90, -45], ccrs.PlateCarree())
ax.coastlines()
gl = ax.gridlines(
draw_labels=True,
lw=0.4,
ylabel_style=dict(path_effects=[PathEffects.withStroke(linewidth=2, foreground="w")]),
xlocs=range(-180,180,45), ylocs=range(-90,0,10),
y_inline=True,
)
cbar = fig.colorbar(im, ax=ax, shrink=0.4, label='Sea Ice Concentration (%)')
ax.set_boundary(mpl.path.Path.circle((0, 0), 3950000*1.2), transform=data_proj)
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论