2

I'm trying to do Azimuthal Equidistant in CF-compliant NetCDF, but I didn't get the correct projection. If you compare the first two information (gdalinfo) you see there is no AUTHORITY and also PROJCS, GEOGCS, DATUM, SPHEROID are unknown.

I followed this link to set the set the projection: Mapping from CF Grid Mapping Attributes to CRS WKT Elements

This is the original projection which I want to save it to NetCDF.

Size is 8000, 8000
Coordinate System is:
PROJCS["Azimuthal_Equidistant",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Azimuthal_Equidistant"],
    PARAMETER["latitude_of_center",53],
    PARAMETER["longitude_of_center",24],
    PARAMETER["false_easting",5837287.81977],
    PARAMETER["false_northing",2121415.69617],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (4800000.000000000000000,1800000.000000000000000)
Pixel Size = (75.000000000000000,-75.000000000000000)

and this is what I got:

Size is 8000, 8000
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["Spheroid",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Azimuthal_Equidistant"],
    PARAMETER["latitude_of_center",53],
    PARAMETER["longitude_of_center",24],
    PARAMETER["false_easting",5837287.81977],
    PARAMETER["false_northing",2121415.69617]]
Origin = (4799602.500000000000000,1800037.500000000000000)
Pixel Size = (75.000000000000000,-75.000000000000000)

here is the code that I extracted the projection:

# Extract Projection info from tiff
ds = gdal.Open('sample.tiff')
prj = ds.GetProjection()
srs = osr.SpatialReference(wkt=prj)


# Create container variable for CRS
crso = nco.createVariable('crs', 'i4')
crso.grid_mapping_name = 'azimuthal_equidistant'
crso.projected_crs_name = srs.GetAttrValue('PROJCS')
crso.geographic_crs_name = srs.GetAttrValue('GEOGCS')
crso.horizontal_datum_name = srs.GetAttrValue('DATUM')
crso.latitude_of_projection_origin = srs.GetProjParm('latitude_of_center')
crso.longitude_of_projection_origin = srs.GetProjParm('longitude_of_center')
crso.false_easting = srs.GetProjParm('false_easting')
crso.false_northing = srs.GetProjParm('false_northing')
crso.longitude_of_prime_meridian = 0.0
crso.semi_major_axis = 6378137.0
crso.inverse_flattening = 298.257223563

Does someone know how to get the correct projection?

2

I found the solution (maybe this will help someone), I have to write the Spatial reference into the CRS variable, then netcdf can recognize all projection info which is defined.

crso.spatial_ref = str(srs)

str(srs) gives me:

 crs#spatial_ref=PROJCS["Azimuthal_Equidistant",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Azimuthal_Equidistant"],
    PARAMETER["latitude_of_center",53],
    PARAMETER["longitude_of_center",24],
    PARAMETER["false_easting",5837287.81977],
    PARAMETER["false_northing",2121415.69617],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.