Geographic Information Systems Asked on August 18, 2021
I am working on a PyQgis application which loads WRF model output as a raster to overlay on a map which is in the Pseudo-Mercator projection. The WRF data is projected to Lambert-Conformal and am using PyProj to retrieve the projection metadata and GDAL to transform the original NetCDF file to GeoTiff. Initially I was getting issues with "gdal_translate" not being able to decipher the Proj4 string I produce, so I decided to use the GDAL API directly to debug; this required converting the Proj4 string to an OGC WKT. I use the QProcess class in the application to run the following:
gdalsrsinfo -o wkt ${MY_PROJ_STR}
where "${MY_PROJ_STR}" is:
"+proj=lcc +lat_1=13.342499732971191 +lat_2=47.317501068115234 +lat_0=30.330001831054688 +lon_0=-86.16000366210938 +x_0=0+y_0=0 +a=6370000 +datum=WGS84 +no_defs"
When I run the "gdalsrsinfo" command from the application I get the following stderr output:
ERROR 1: ERROR - failed to load SRS definition from "+proj=lcc +lat_1=13.342499732971191 +lat_2=47.317501068115234 +lat_0=30.330001831054688 +lon_0=-86.16000366210938 +x_0=0+y_0=0 +a=6370000 +datum=WGS84 +no_defs"
yet when testing it from the command line, it get the WKT translation:
PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",13.34249973297119],PARAMETER["standard_parallel_2",47.31750106811523],PARAMETER["latitude_of_origin",30.33000183105469],PARAMETER["central_meridian",-86.16000366210938],PARAMETER["false_easting",0],PARAMETER["false_northing",0]]
I’m not sure if it’s the PROJCS="unnamed" part that’s throwing it off. I checked to make sure the versions of GDAL used in Python and from the command line weren’t in conflict, and they seem to be the same:
print(gdal.VersionInfo())
>> 2030200
gdalinfo --version
]$ GDAL 2.3.2, released 2018/09/21
My other system specs for reference: am using Python 3.7, PyQt5 5.13.2, pyproj 3.0.0.post1, PROJ 5.2.0, and QGIS 5.10. The operating system is Fedora 31. I’m not sure if the syntax of the PROJ string is correct or if I am missing a step in producing the correct reference system. I’ve also added "GDAL_DATA" and "PROJ_LIB" environmental variables to my ~/.bash_profile.
I'm not sure where your issue is as I get the same result using below samples. PS: I suppose you will use Python API to call GDAL and not quick execution of GDAL binary within Python using subprocess
or os.open
With Python and GDAL combined
import osr
srs_def = "+proj=lcc +lat_1=13.342499732971191 +lat_2=47.317501068115234 +lat_0=30.330001831054688 +lon_0=-86.16000366210938 +x_0=0+y_0=0 +a=6370000 +datum=WGS84 +no_defs"
proj = osr.SpatialReference()
proj.ImportFromProj4(srs_def)
print(proj.ExportToWkt())
print(proj.ExportToPrettyWkt())
With command line
gdalsrsinfo "+proj=lcc +lat_1=13.342499732971191 +lat_2=47.317501068115234 +lat_0=30.330001831054688 +lon_0=-86.16000366210938 +x_0=0+y_0=0 +a=6370000 +datum=WGS84 +no_defs"
Correct answer by ThomasG77 on August 18, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP