如何在Python中绘制netCDF海冰范围数据集？

huangapple
go评论131阅读模式
英文:

How do I plot netCDF sea ice extent dataset in Python?

问题

我计划使用Python绘制海洋冰盖在南极附近的分布图，使用NOAA-NSIDC被动微波海冰浓度第4版本数据集。我正在使用以下代码：

  1. import xarray as xr
  2. import numpy as np
  3. import cmocean
  4. import cartopy.crs as ccrs
  5. import matplotlib.pyplot as plt
  6. import cartopy.feature as cfeature
  7. from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
  9. # 读取netCDF文件
  10. ds = xr.open_dataset('seaice_conc_monthly_sh_202212_f17_v04r00.nc')
  12. # 提取第一个时间步的海冰范围数据
  13. sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values
  15. # 提取x和y网格坐标
  16. xgrid = ds.xgrid.values
  17. ygrid = ds.ygrid.values
  19. # 创建图形
  20. fig = plt.figure(figsize=(10, 10))
  21. ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
  22. ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
  23. im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=ccrs.PlateCarree(), cmap='cmo.ice')
  24. ax.gridlines()
  25. ax.coastlines()
  26. gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
  27. gl.top_labels = False
  28. gl.right_labels = False
  29. gl.xformatter = LONGITUDE_FORMATTER
  30. gl.yformatter = LATITUDE_FORMATTER
  31. cbar = fig.colorbar(im, ax=ax, shrink=0.4)
  32. cbar.ax.set_ylabel('Sea Ice Concentration (%)', fontsize=16)
  33. plt.show()

然而，我遇到了以下错误：

  1. 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

我想要获得类似于此图的绘图：
如何在Python中绘制netCDF海冰范围数据集？

英文:

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:

  1. >>> import xarray as xr
  2. >>> import numpy as np
  3. >>> import cmocean
  4. >>> import cartopy.crs as ccrs
  5. >>> import matplotlib.pyplot as plt
  6. >>> import cartopy.feature as cfeature
  7. >>> from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
  8. >>> # Read in the netCDF file
  9. >>>
  10. >>> ds = xr.open_dataset('seaice_conc_monthly_sh_202212_f17_v04r00.nc')
  11. >>> # Extract the sea ice extent data for the first time step
  12. >>>
  13. >>> sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values
  14. >>> # Extract the x and y grid coordinates
  15. >>>
  16. >>> xgrid = ds.xgrid.values
  17. >>> ygrid = ds.ygrid.values
  18. >>> # Create the plot
  19. >>>
  20. >>> fig = plt.figure(figsize=(10, 10))
  21. >>> ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
  22. >>> ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
  23. >>> im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=ccrs.PlateCarree(), cmap='cmo.ice')
  24. Traceback (most recent call last):
  25. File "<stdin>", line 1, in <module>
  26. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 318, in wrapper
  27. return func(self, *args, **kwargs)
  28. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 1787, in pcolormesh
  29. result = self._wrap_quadmesh(result, **kwargs)
  30. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 1970, in _wrap_quadmesh
  31. pcolor_col = self.pcolor(coords[..., 0], coords[..., 1],
  32. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 318, in wrapper
  33. return func(self, *args, **kwargs)
  34. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 2015, in pcolor
  35. limits = result.get_datalim(self.transData)
  36. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 266, in get_datalim
  37. paths = [transform.transform_path_non_affine(p) for p in paths]
  38. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 266, in <listcomp>
  39. paths = [transform.transform_path_non_affine(p) for p in paths]
  40. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\transforms.py", line 2431, in transform_path_non_affine
  41. return self._a.transform_path_non_affine(path)
  42. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 190, in transform_path_non_affine
  43. proj_geom = self.target_projection.project_geometry(
  44. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 808, in project_geometry
  45. return getattr(self, method_name)(geometry, src_crs)
  46. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 963, in _project_polygon
  47. return self._rings_to_multi_polygon(rings, is_ccw)
  48. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 1224, in _rings_to_multi_polygon
  49. multi_poly = sgeom.MultiPolygon(polygon_bits)
  50. File "C:\Users\ashif\miniconda3\lib\site-packages\shapely\geometry\multipolygon.py", line 79, in __new__
  51. shell = ob[0]
  52. TypeError: 'GeometryCollection' object is not subscriptable
  53. >>> ax.gridlines()
  54. <cartopy.mpl.gridliner.Gridliner object at 0x000001EE3C35B8E0>
  55. >>> ax.coastlines()
  56. <cartopy.mpl.feature_artist.FeatureArtist object at 0x000001EE3C31FE80>
  57. >>> gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
  58. >>> gl.top_labels = False
  59. >>> gl.right_labels = False
  60. >>> gl.xformatter = LONGITUDE_FORMATTER
  61. >>> gl.yformatter = LATITUDE_FORMATTER
  62. >>> cbar = fig.colorbar(im, ax=ax, shrink=0.4)
  63. Traceback (most recent call last):
  64. File "<stdin>", line 1, in <module>
  65. NameError: name 'im' is not defined. Did you mean: 'id'?
  66. >>> cbar.ax.set_ylabel('Sea Ice Concentration (%)', fontsize=16)
  67. Traceback (most recent call last):
  68. File "<stdin>", line 1, in <module>
  69. NameError: name 'cbar' is not defined
  70. >>> plt.show()
  71. Traceback (most recent call last):
  72. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\backends\backend_qt.py", line 468, in _draw_idle
  73. self.draw()
  74. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\backends\backend_agg.py", line 400, in draw
  75. self.figure.draw(self.renderer)
  76. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 95, in draw_wrapper
  77. result = draw(artist, renderer, *args, **kwargs)
  78. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
  79. return draw(artist, renderer)
  80. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\figure.py", line 3140, in draw
  81. mimage._draw_list_compositing_images(
  82. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\image.py", line 131, in _draw_list_compositing_images
  83. a.draw(renderer)
  84. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
  85. return draw(artist, renderer)
  86. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 538, in draw
  87. return super().draw(renderer=renderer, **kwargs)
  88. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
  89. return draw(artist, renderer)
  90. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\axes\_base.py", line 3064, in draw
  91. mimage._draw_list_compositing_images(
  92. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\image.py", line 131, in _draw_list_compositing_images
  93. a.draw(renderer)
  94. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
  95. return draw(artist, renderer)
  96. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 972, in draw
  97. super().draw(renderer)
  98. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\artist.py", line 72, in draw_wrapper
  99. return draw(artist, renderer)
  100. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 351, in draw
  101. transform, offset_trf, offsets, paths = self._prepare_points()
  102. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 328, in _prepare_points
  103. paths = [transform.transform_path_non_affine(path)
  104. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\collections.py", line 328, in <listcomp>
  105. paths = [transform.transform_path_non_affine(path)
  106. File "C:\Users\ashif\miniconda3\lib\site-packages\matplotlib\transforms.py", line 2431, in transform_path_non_affine
  107. return self._a.transform_path_non_affine(path)
  108. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\mpl\geoaxes.py", line 190, in transform_path_non_affine
  109. proj_geom = self.target_projection.project_geometry(
  110. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 808, in project_geometry
  111. return getattr(self, method_name)(geometry, src_crs)
  112. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 963, in _project_polygon
  113. return self._rings_to_multi_polygon(rings, is_ccw)
  114. File "C:\Users\ashif\miniconda3\lib\site-packages\cartopy\crs.py", line 1224, in _rings_to_multi_polygon
  115. multi_poly = sgeom.MultiPolygon(polygon_bits)
  116. File "C:\Users\ashif\miniconda3\lib\site-packages\shapely\geometry\multipolygon.py", line 79, in __new__
  117. shell = ob[0]
  118. TypeError: 'GeometryCollection' object is not subscriptable

However I am facing the following error:

  1. 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:

  1. Attributes:
  2. long_name:            NOAA/NSIDC Climate Data Record of Passive Microwave...
  3. standard_name:        sea_ice_area_fraction
  4. units:                1
  5. flag_values:          [251 252 253 254 255]
  6. flag_meanings:        pole_hole lakes coastal land_mask missing_data
  7. datum:                +ellps=urn:ogc:def:crs:EPSG::4326
  8. grid_mapping:         projection
  9. reference:            https://nsidc.org/data/g02202/versions/4/
  10. ancillary_variables:  stdev_of_cdr_seaice_conc_monthly qa_of_cdr_seaice_c...
  11. valid_range:          [  0 100], 'nsidc_bt_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
  12. [104912 values with dtype=float32]
  13. Attributes:
  14. long_name:      Passive Microwave Monthly Southern Hemisphere Sea Ice Con...
  15. standard_name:  sea_ice_area_fraction
  16. units:          1
  17. flag_values:    [251 252 253 254 255]
  18. flag_meanings:  pole_hole unused coastal land_mask missing_data
  19. datum:          +ellps=urn:ogc:def:crs:EPSG::4326
  20. grid_mapping:   projection
  21. valid_range:    [  0 100], 'nsidc_nt_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
  22. [104912 values with dtype=float32]
  23. Attributes:
  24. long_name:      Passive Microwave Monthly Southern Hemisphere Sea Ice Con...
  25. standard_name:  sea_ice_area_fraction
  26. units:          1
  27. flag_values:    [251 252 253 254 255]
  28. flag_meanings:  pole_hole unused coastal land_mask missing_data
  29. datum:          +ellps=urn:ogc:def:crs:EPSG::4326
  30. grid_mapping:   projection
  31. valid_range:    [  0 100], 'projection': <xarray.Variable ()>
  32. [1 values with dtype=|S1]
  33. Attributes: (12/23)
  34. grid_boundary_top_projected_y:          4350000.0
  35. grid_boundary_bottom_projected_y:       -3950000.0
  36. grid_boundary_right_projected_x:        3950000.0
  37. grid_boundary_left_projected_x:         -3950000.0
  38. parent_grid_cell_row_subset_start:      0.0
  39. parent_grid_cell_row_subset_end:        332.0
  40. ...                                     ...
  41. scaling_factor:                         1.0
  42. false_easting:                          0.0
  43. false_northing:                         0.0
  44. semimajor_radius:                       6378273.0
  45. semiminor_radius:                       6356889.449
  46. units:                                  meters, 'qa_of_cdr_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
  47. [104912 values with dtype=float32]
  48. Attributes:
  49. long_name:      Passive Microwave Monthly Southern Hemisphere Sea Ice Con...
  50. standard_name:  sea_ice_area_fraction status_flag
  51. flag_meanings:  Average_concentration_exceeds_0.15 Average_concentration_...
  52. datum:          +ellps=urn:ogc:def:crs:EPSG::4326
  53. grid_mapping:   projection
  54. flag_masks:     [  1   2   4   8  32  64 128]
  55. valid_range:    [  1 255], 'stdev_of_cdr_seaice_conc_monthly': <xarray.Variable (tdim: 1, y: 332, x: 316)>
  56. [104912 values with dtype=float32]
  57. Attributes:
  58. long_name:     Passive Microwave Monthly Southern Hemisphere Sea Ice Conc...
  59. datum:         +ellps=urn:ogc:def:crs:EPSG::4326
  60. grid_mapping:  projection
  61. valid_range:   [0. 1.], 'time': <xarray.Variable (tdim: 1)>
  62. [1 values with dtype=datetime64[ns]]
  63. Attributes:
  64. standard_name:  time
  65. long_name:      ANSI date
  66. axis:           T, 'xgrid': <xarray.Variable (x: 316)>
  67. array([-3937500., -3912500., -3887500., ...,  3887500.,  3912500.,  3937500.],
  68. dtype=float32)
  69. Attributes:
  70. valid_range:    [-3950000.  3950000.]
  71. units:          meters
  72. long_name:      projection_grid_x_centers
  73. standard_name:  projection_x_coordinate
  74. axis:           X, 'ygrid': <xarray.Variable (y: 332)>
  75. array([ 4337500.,  4312500.,  4287500., ..., -3887500., -3912500., -3937500.],
  76. dtype=float32)
  77. Attributes:
  78. valid_range:    [-3950000.  4350000.]
  79. units:          meters
  80. long_name:      projection_grid_y_centers
  81. standard_name:  projection_y_coordinate
  82. axis:           Y}))

The following variables are in the dataset:

  1. cdr_seaice_conc_monthly           (tdim, y, x) float32 ...
  2. nsidc_bt_seaice_conc_monthly      (tdim, y, x) float32 ...
  3. nsidc_nt_seaice_conc_monthly      (tdim, y, x) float32 ...
  4. projection                        |S1 ...
  5. qa_of_cdr_seaice_conc_monthly     (tdim, y, x) float32 ...
  6. stdev_of_cdr_seaice_conc_monthly  (tdim, y, x) float32 ...

I would like to get a plot similar to this:

如何在Python中绘制netCDF海冰范围数据集？

答案1

得分: 0

主要问题在于您指定了数据的投影为ccrs.PlateCarree()，这显然是不正确的。

也许有比下面示例中解析更好的方法，但它捕捉到了要点：

  1. import xarray as xr
  2. import numpy as np
  3. import cartopy.crs as ccrs
  4. import matplotlib.pyplot as plt
  5. import matplotlib as mpl
  6. import matplotlib.patheffects as PathEffects
  8. with xr.open_dataset("seaice_conc_monthly_sh_202212_f17_v04r00.nc") as ds:
  10.     sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values * 100
  12.     data_proj = ccrs.epsg(ds.projection.attrs["srid"].split("::")[-1])
  13.     
  14.     flags = ds.cdr_seaice_conc_monthly.attrs["flag_values"]
  15.     
  16.     invalid = np.isin(sea_ice_extent, flags)
  17.     sea_ice_extent[invalid] = np.nan
  18.     
  19.     xgrid = ds.xgrid.values
  20.     ygrid = ds.ygrid.values
  22. cmap = mpl.colormaps["BuPu"].copy()
  23. cmap.set_bad("grey", alpha=.2)
  25. fig, ax = plt.subplots(
  26.     figsize=(8,8), facecolor="w",
  27.     subplot_kw=dict(projection=ccrs.SouthPolarStereo()),
  28. )
  30. im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=data_proj, cmap=cmap)
  32. ax.set_extent([-180, 180, -90, -45], ccrs.PlateCarree())
  33. ax.coastlines()
  35. gl = ax.gridlines(
  36.     draw_labels=True, 
  37.     lw=0.4, 
  38.     ylabel_style=dict(path_effects=[PathEffects.withStroke(linewidth=2, foreground="w")]),
  39.     xlocs=range(-180,180,45), ylocs=range(-90,0,10),
  40.     y_inline=True,
  41. )
  42. cbar = fig.colorbar(im, ax=ax, shrink=0.4, label='Sea Ice Concentration (%)')
  44. ax.set_boundary(mpl.path.Path.circle((0, 0), 3950000*1.2), transform=data_proj)

如何在Python中绘制netCDF海冰范围数据集？

英文:

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:

  1. import xarray as xr
  2. import numpy as np
  3. import cartopy.crs as ccrs
  4. import matplotlib.pyplot as plt
  5. import matplotlib as mpl
  6. import matplotlib.patheffects as PathEffects
  8. with xr.open_dataset("seaice_conc_monthly_sh_202212_f17_v04r00.nc") as ds:
  10.     sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values * 100
  12.     data_proj = ccrs.epsg(ds.projection.attrs["srid"].split("::")[-1])
  13.     
  14.     flags = ds.cdr_seaice_conc_monthly.attrs["flag_values"]
  15.     
  16.     invalid = np.isin(sea_ice_extent, flags)
  17.     sea_ice_extent[invalid] = np.nan
  18.     
  19.     xgrid = ds.xgrid.values
  20.     ygrid = ds.ygrid.values
  22. cmap = mpl.colormaps["BuPu"].copy()
  23. cmap.set_bad("grey", alpha=.2)
  25. fig, ax = plt.subplots(
  26.     figsize=(8,8), facecolor="w",
  27.     subplot_kw=dict(projection=ccrs.SouthPolarStereo()),
  28. )
  30. im = ax.pcolormesh(xgrid, ygrid, sea_ice_extent, transform=data_proj, cmap=cmap)
  32. ax.set_extent([-180, 180, -90, -45], ccrs.PlateCarree())
  33. ax.coastlines()
  35. gl = ax.gridlines(
  36.     draw_labels=True, 
  37.     lw=0.4, 
  38.     ylabel_style=dict(path_effects=[PathEffects.withStroke(linewidth=2, foreground="w")]),
  39.     xlocs=range(-180,180,45), ylocs=range(-90,0,10),
  40.     y_inline=True,
  41. )
  42. cbar = fig.colorbar(im, ax=ax, shrink=0.4, label='Sea Ice Concentration (%)')
  44. ax.set_boundary(mpl.path.Path.circle((0, 0), 3950000*1.2), transform=data_proj)

如何在Python中绘制netCDF海冰范围数据集？

huangapple
  • 本文由 发表于 2023年5月31日 23:24:58
  • 转载请务必保留本文链接：https://go.coder-hub.com/76375088.html
  • cartopy
  • netcdf
  • python
  • python-xarray
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定