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
<https://github.com/cf-convention/cf-conventions/wiki/Mapping-from-CF-Grid-Mapping-Attributes-to-CRS-WKT-Elements>
This is the original projection which I want to save it to NetCDF.
|
Sizeis8000,8000
CoordinateSystemis:
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)
PixelSize=(75.000000000000000,-75.000000000000000)
|
and this is what I got:
|
Sizeis8000,8000
CoordinateSystemis:
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)
PixelSize=(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?
