[gdal-dev] GetProjection returns empty string, python 3.7

Nicolas Cadieux njacadieux.gitlab at gmail.com
Mon Oct 26 10:53:39 PDT 2020


Hi,

In case mail question was not clear, I have posted the questions about 
this on stackexchange.  I posted a full code that you can run exposing 
my problem and question.

Thanks

Nicolas

https://gis.stackexchange.com/questions/377567/cannot-get-projection-from-raster-dataset-using-getprojection

On 2020-10-25 9:31 p.m., Nicolas Cadieux wrote:
> Hi,
>
> With the following code, I get an empty string indicating the 
> projection is not valid.
>
> from osgeo import gdal, osr
> raster_ds = gdal.Open(r"C:\temp\180922_WTE3.tif")
> target_ds = 
> driver.Create(r"c:\temp\output.tif",xsize=raster_ds_ncol,ysize=raster_ds_nrow,bands 
> = 1,eType = gdal.GDT_Float32)
> (array is a numpty array.)
> raster_srs = osr.SpatialReference()
> raster_srs.ImportFromWkt(raster_ds.GetProjectionRef())
> target_ds.SetProjection(raster_srs.ExportToWkt())
> target_ds.GetRasterBand(1).WriteArray(array)
> raster_ds = None #close dataset
> target_ds=None
>
> Below is the result of gdal info from qgis. File appears to have a 
> valid projection and is properly georeferenced in QGIS using other 
> data sets.   Is this projection wrong or am I missing something in my 
> python code?  The goal is to extract the projection from raster_ds set 
> in order to apply to target_ds.  I can create the file, apply a 
> geotransform but can't get the projection.  Can you help?
>
> Thanks/merci!
>
> Nicolas
>
> QGIS version: 3.14.16-Pi
> QGIS code revision: df27394552
> Qt version: 5.11.2
> GDAL version: 3.0.4
> GEOS version: 3.8.1-CAPI-1.13.3
> PROJ version: Rel. 6.3.2, May 1st, 2020
> Processing algorithm…
> Algorithm 'Raster information' starting…
> Input parameters:
> { 'EXTRA' : '', 'INPUT' : 'C:/temp/180922_WTE3.tif', 'MIN_MAX' : 
> False, 'NOGCP' : False, 'NO_METADATA' : False, 'OUTPUT' : 
> 'TEMPORARY_OUTPUT', 'STATS' : False }
>
> GDAL command:
> gdalinfo C:/temp/180922_WTE3.tif
> GDAL command output:
> Warning 1: TIFFReadDirectory:Sum of Photometric type-related color 
> channels and ExtraSamples doesn't match SamplesPerPixel. Defining 
> non-color channels as ExtraSamples.
>
> Driver: GTiff/GeoTIFF
> Files: C:/temp/180922_WTE3.tif
> Size is 1194, 2783
> Coordinate System is:
> PROJCRS["WGS 84 / UTM zone 18N",
> BASEGEOGCRS["WGS 84",
> DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433]],
> ID["EPSG",4326]],
> CONVERSION["UTM zone 18N",
> METHOD["Transverse Mercator",
> ID["EPSG",9807]],
> PARAMETER["Latitude of natural origin",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8801]],
> PARAMETER["Longitude of natural origin",-75,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8802]],
> PARAMETER["Scale factor at natural origin",0.9996,
> SCALEUNIT["unity",1],
> ID["EPSG",8805]],
> PARAMETER["False easting",500000,
> LENGTHUNIT["metre",1],
> ID["EPSG",8806]],
> PARAMETER["False northing",0,
> LENGTHUNIT["metre",1],
> ID["EPSG",8807]]],
> CS[Cartesian,2],
> AXIS["(E)",east,
> ORDER[1],
> LENGTHUNIT["metre",1]],
> AXIS["(N)",north,
> ORDER[2],
> LENGTHUNIT["metre",1]],
> USAGE[
> SCOPE["unknown"],
> AREA["World - N hemisphere - 78°W to 72°W - by country"],
> BBOX[0,-78,84,-72]],
> ID["EPSG",32618]]
> Data axis to CRS axis mapping: 1,2
> Origin = (459417.200000000011642,5028584.700000000186265)
> Pixel Size = (0.050000000000000,-0.050000000000000)
> Metadata:
> AREA_OR_POINT=Area
> TIFFTAG_XRESOLUTION=1
> TIFFTAG_YRESOLUTION=1
> Image Structure Metadata:
> INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left ( 459417.200, 5028584.700) ( 75d31' 7.03"W, 45d24'34.58"N)
> Lower Left ( 459417.200, 5028445.550) ( 75d31' 6.99"W, 45d24'30.07"N)
> Upper Right ( 459476.900, 5028584.700) ( 75d31' 4.28"W, 45d24'34.59"N)
> Lower Right ( 459476.900, 5028445.550) ( 75d31' 4.24"W, 45d24'30.08"N)
> Center ( 459447.050, 5028515.125) ( 75d31' 5.63"W, 45d24'32.33"N)
> Band 1 Block=1194x1 Type=UInt16, ColorInterp=Red
> ...
> Band 288 Block=1194x1 Type=UInt16, ColorInterp=Undefined
>
>
>


More information about the gdal-dev mailing list