Use the proj.4 library and pyproj
On Feb 10, 2016 5:23 AM, "icloud" <khakpour.m@xxxxxx> wrote:
>
> 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.
>
> 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?
>
>
>
>
> _______________________________________________
> netcdfgroup mailing list
> netcdfgroup@xxxxxxxxxxxxxxxx
> For list information or to unsubscribe, visit:
> http://www.unidata.ucar.edu/mailing_lists/
>
