Hi,

Looks like it's an environment thing... Both versions of the code work with 3.1.4.

Thanks for your help!

Nicolas

On 2020-10-26 6:31 p.m., Paul Harwood wrote:
Me too.

But which conda package for gdal are you using - there are a lot of ones out there.

Are you using conda install -c conda-forge gdal ?

It should be version 3.1.4 at the moment ?

Also - I don't know if this is part of it but ...

"srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("NAD27")
dst_ds.SetProjection(srs.ExportToWkt())
"

why not just

"srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("NAD27")
dst_ds.SetSpatialRef(srs)
"
??

On Mon, 26 Oct 2020 at 22:11, Nicolas Cadieux <njacadieux.git...@gmail.com <mailto:njacadieux.git...@gmail.com>> wrote:

    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
    <njacadieux.git...@gmail.com
    <mailto:njacadieux.git...@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.git...@gmail.com
        <mailto:njacadieux.git...@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@lists.osgeo.org <mailto:gdal-dev@lists.osgeo.org>
            https://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to