时间戳

首页 > 解决方案 > 从 proj4 参数构造 cartopy CRS

python - 从 proj4 参数构造 cartopy CRS

问题描述

我下载了 MODIS .hdf 文件。用 xarray 加载它，它给了我一个属性ds.Proj4String == ' +a=6378137.0 +b=6356752.3142451793 +no_defs +proj=latlong\n'。

如何使用该字符串将原始坐标转换为例如纬度/经度？

这是数据的样子：

<xarray.Dataset>
Dimensions:                  (XDim:mod06: 1503, YDim:mod06: 833)
Dimensions without coordinates: XDim:mod06, YDim:mod06
Data variables:
    Cloud_Optical_Thickness  (YDim:mod06, XDim:mod06) float32 ...
Attributes:
    HDFEOSVersion:     HDFEOS_V2.19
    StructMetadata.0:  GROUP=SwathStructure\nEND_GROUP=SwathStructure\nGROUP=...
    CoreMetadata:      \nGROUP = INVENTORYMETADATA\n  GROUPTYPE = MASTERGROUP...
    ArchiveMetadata:   GROUP = ARCHIVEDMETADATA\n  GROUPTYPE = MASTERGROUP\n\...
    Proj4String:        +a=6378137.0 +b=6356752.3142451793 +no_defs +proj=lat...

文档proj4_params告诉我作为键值对传递，所以我去了

ss = ds.Proj4String.split()
proj = {}
for s in ss:
    k = s.split('=')
    if len(k)==2:
        proj[k[0][1:]] = k[1]
print(proj)

并到达{'a': '6378137.0', 'b': '6356752.3142451793', 'proj': 'latlong'}，但ccrs.CRS(proj)仍会引发有关globe缺少属性的错误。

标签： pythongiscartopyproj

解决方案

尚不支持仅从 proj.4 字符串创建 CartoPy 投影，请参阅此 PR。您使用的 CRS 类本身并不足以在 CartoPy 中创建功能齐全的投影，并且实际上被设计为其他投影的基类。

就您得到的错误而言，globe无论传入什么，它都是必需的参数proj4_params。您可以使用以下内容创建Globe实例：

globe = ccrs.Globe(semimajor_axis=proj.pop('a'), semiminor_axis=proj.pop('b'))

推荐阅读