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

Paul Harwood runette at gmail.com
Mon Oct 26 14:04:04 PDT 2020 

GetSpatialRef is OSR based and it is certainly working for me against GDAL
Raster datasets - it would not be a method against the GDAL Dataset if it
was only an OGR thing - OGR has its own objects.

This code is working for me at this very minute and giving me the name of
the CRS for a gtiff

dataset = gdal.Open(str(filepath))
        if dataset:
            proj = "None"
            crs = dataset.GetSpatialRef()
            if crs:
                proj = crs.GetName()
            return {"type": "gdal", "driver":
dataset.GetDriver().ShortName, "proj": proj}

On Mon, 26 Oct 2020 at 20:08, Nicolas Cadieux <njacadieux.gitlab at gmail.com>
wrote:

> Hi,
>
> Thanks for helping.  I will put my code at the end.  The problem seems to
> be between gdal python lib 2.2.2 and 3.0.2.  Under 2.2.2 and < 3,
> GetProjection() works well but return nothing with gdal python lib 3.0.2.
> It's my understanding that .GetSpatialRef is for OGR data (vector).  Am I
> wrong?
>
> I was using GetProjection() because that is the method used in the
> raster_api tutorial.  (https://gdal.org/tutorials/raster_api_tut.html)
>
> I understand that gdal 3 now relies on proj 6 for the SpatialRef.  I'am
> trying to figure out if is a problem with proj or with gdal...
>
> Nicolas
>
>
> code:
>
> # -*- coding: utf-8 -*-
> """
> Created on Sun Oct 25 15:09:05 2020
>
> @author: Nicolas
> """
> #
> https://stackoverflow.com/questions/37648439/simplest-way-to-save-array-into-raster-file-in-python
>
>
> from osgeo import gdal, osr, ogr
> import numpy
>
> print (gdal.__version__)
>
> # https://gdal.org/tutorials/raster_api_tut.html
> fileformat = "GTiff"
> driver = gdal.GetDriverByName(fileformat)
> metadata = driver.GetMetadata()
> if metadata.get(gdal.DCAP_CREATE) == "YES":
>     print("Driver {} supports Create() method.".format(fileformat))
>
> if metadata.get(gdal.DCAP_CREATECOPY) == "YES":
>     print("Driver {} supports CreateCopy() method.".format(fileformat))
> dst_ds = driver.Create(r"c:\temp\gdal.tif", xsize=512, ysize=512,
>                     bands=1, eType=gdal.GDT_Byte)
> dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
> srs = osr.SpatialReference()
> srs.SetUTM(11, 1)
> srs.SetWellKnownGeogCS("NAD27")
> dst_ds.SetProjection(srs.ExportToWkt())
> print ('srs = ',srs)# this is good
> raster = numpy.zeros((512, 512), dtype=numpy.uint8)
> dst_ds.GetRasterBand(1).WriteArray(raster)
> # Once we're done, close properly the dataset
> srs = None
> dst_ds = None #srs is file and well georeferenced in Qgis.
>
> #
> ---------------------------------------------------------------------------
> # open dataset try to read srs
> #
> ---------------------------------------------------------------------------
> raster_ds = gdal.Open(r"C:\temp\gdal.tif")
> # first try
> print("Projection is {}".format(raster_ds.GetProjection()))
>
> # second try
> srs = osr.SpatialReference()
> srs.ImportFromWkt(raster_ds.GetProjectionRef())
> print ('srs =', srs)
> # thrid try
> srs = osr.SpatialReference()
> srs = raster_ds.GetProjection()
> print('srs =', srs)
> srs = raster_ds.GetProjectionRef()
> # forth try
> # from gdal_info
> https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/gdalinfo.py
> pszProjection = raster_ds.GetProjectionRef()
> print(pszProjection)
> if pszProjection is not None:
>     hSRS = osr.SpatialReference()
>     if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None:
>         pszPrettyWkt = hSRS.ExportToPrettyWkt(False)
>         print("Coordinate System is:\n%s" % pszPrettyWkt)
>     else:
>         print("Coordinate System is `%s'" % pszProjection)
>
>
> On 2020-10-26 3:44 p.m., Paul Harwood wrote:
>
> I am not sure why you are using GetProjection and not GetSpatialRef ?
> Although, the code snippet you included in your email does not seem to use
> either?
>
> I certainly was in the process of working on some code that uses
> .GeoSptialRef very successfully on GDAL3.1.3 - so I am sure that it is
> working!
>
> On Mon, 26 Oct 2020 at 18:36, Nicolas Cadieux <njacadieux.gitlab at gmail.com>
> wrote:
>
>> New info!
>>
>> GetProjection works with gdal python 2.2.2 and 2.3.3 but not with with
>> gdal 3.0.2.  Is this a change in the gdal library?
>>
>> Nicolas
>>
>> On 2020-10-26 1:53 p.m., Nicolas Cadieux wrote:
>> > 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
>> >>
>> >>
>> >>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201026/9c57c273/attachment.html>

More information about the gdal-dev mailing list