CODEWITHSUNDEEP
pythongiscartopyproj
doppler
doppler

Reputation: 997

Construct cartopy CRS from proj4 parameters

I downloaded a MODIS .hdf file. Loading it with xarray, it gives me an attribute ds.Proj4String == ' +a=6378137.0 +b=6356752.3142451793 +no_defs +proj=latlong\n'.

How can I use that string to convert the raw coordinates to e.g. lat/lon?

This is what the data look like:

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

The docs tell me to pass proj4_params as key-value pairs, so I go

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

and arrive at {'a': '6378137.0', 'b': '6356752.3142451793', 'proj': 'latlong'}, but ccrs.CRS(proj) still throws an error about attribute globe missing.

Upvotes: 3

Views: 1725

Answers (1)

DopplerShift
DopplerShift

Reputation: 5843

Creating CartoPy projections solely from proj.4 strings is not supported yet, see this PR. The CRS class you're using really is insufficient by itself to create a fully functioning projection in CartoPy, and is really designed to be a base class for the other projections.

As far as the error you're getting, globe is a required parameter regardless of what's being passed in proj4_params. You could create the Globe instance with something like:

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

Upvotes: 6

Related Questions