Hi,

Comes back as None for me... :(  I'am working with Anaconda. Could it be my environment?

Nicolas

On 2020-10-26 5:04 p.m., Paul Harwood wrote:
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 <[email protected] <mailto:[email protected]>> 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
    <[email protected]
    <mailto:[email protected]>> 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
        [email protected] <mailto:[email protected]>
        https://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to