<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
  </head>
  <body>
    <p>Hi,</p>
    <p>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?  <br>
    </p>
    <p>I was using GetProjection() because that is the method used in
      the raster_api tutorial. 
      (<a class="moz-txt-link-freetext" href="https://gdal.org/tutorials/raster_api_tut.html">https://gdal.org/tutorials/raster_api_tut.html</a>)</p>
    <p>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...<br>
    </p>
    <p>Nicolas<br>
    </p>
    <p><br>
    </p>
    <p>code:</p>
    <p># -*- coding: utf-8 -*-<br>
      """<br>
      Created on Sun Oct 25 15:09:05 2020<br>
      <br>
      @author: Nicolas<br>
      """<br>
#<a class="moz-txt-link-freetext" href="https://stackoverflow.com/questions/37648439/simplest-way-to-save-array-into-raster-file-in-python">https://stackoverflow.com/questions/37648439/simplest-way-to-save-array-into-raster-file-in-python</a><br>
      <br>
      <br>
      from osgeo import gdal, osr, ogr<br>
      import numpy<br>
      <br>
      print (gdal.__version__)<br>
       <br>
      # <a class="moz-txt-link-freetext" href="https://gdal.org/tutorials/raster_api_tut.html">https://gdal.org/tutorials/raster_api_tut.html</a><br>
      fileformat = "GTiff"<br>
      driver = gdal.GetDriverByName(fileformat)<br>
      metadata = driver.GetMetadata()<br>
      if metadata.get(gdal.DCAP_CREATE) == "YES":<br>
          print("Driver {} supports Create()
      method.".format(fileformat))<br>
      <br>
      if metadata.get(gdal.DCAP_CREATECOPY) == "YES":<br>
          print("Driver {} supports CreateCopy()
      method.".format(fileformat))<br>
      dst_ds = driver.Create(r"c:\temp\gdal.tif", xsize=512, ysize=512,<br>
                          bands=1, eType=gdal.GDT_Byte)<br>
      dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])<br>
      srs = osr.SpatialReference()<br>
      srs.SetUTM(11, 1)<br>
      srs.SetWellKnownGeogCS("NAD27")<br>
      dst_ds.SetProjection(srs.ExportToWkt())<br>
      print ('srs = ',srs)# this is good<br>
      raster = numpy.zeros((512, 512), dtype=numpy.uint8)<br>
      dst_ds.GetRasterBand(1).WriteArray(raster)<br>
      # Once we're done, close properly the dataset<br>
      srs = None<br>
      dst_ds = None #srs is file and well georeferenced in Qgis.<br>
      <br>
      #
---------------------------------------------------------------------------<br>
      # open dataset try to read srs<br>
      #
---------------------------------------------------------------------------<br>
      raster_ds = gdal.Open(r"C:\temp\gdal.tif")<br>
      # first try<br>
      print("Projection is {}".format(raster_ds.GetProjection()))<br>
      <br>
      # second try<br>
      srs = osr.SpatialReference()<br>
      srs.ImportFromWkt(raster_ds.GetProjectionRef())<br>
      print ('srs =', srs)<br>
      # thrid try<br>
      srs = osr.SpatialReference()<br>
      srs = raster_ds.GetProjection()<br>
      print('srs =', srs)<br>
      srs = raster_ds.GetProjectionRef()<br>
      # forth try<br>
      # from gdal_info
<a class="moz-txt-link-freetext" href="https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/gdalinfo.py">https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/gdalinfo.py</a><br>
      pszProjection = raster_ds.GetProjectionRef()<br>
      print(pszProjection)<br>
      if pszProjection is not None:<br>
          hSRS = osr.SpatialReference()<br>
          if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None:<br>
              pszPrettyWkt = hSRS.ExportToPrettyWkt(False)<br>
              print("Coordinate System is:\n%s" % pszPrettyWkt)<br>
          else:<br>
              print("Coordinate System is `%s'" % pszProjection)<br>
      <br>
    </p>
    <p><br>
    </p>
    <div class="moz-cite-prefix">On 2020-10-26 3:44 p.m., Paul Harwood
      wrote:<br>
    </div>
    <blockquote type="cite"
cite="mid:CAE8nN5POG7WXyPyLxazxL2FXEnb=669DjVW2LjB_ATeNFME4ww@mail.gmail.com">
      <meta http-equiv="content-type" content="text/html; charset=UTF-8">
      <div dir="ltr">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?<br>
        <br>
        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!</div>
      <br>
      <div class="gmail_quote">
        <div dir="ltr" class="gmail_attr">On Mon, 26 Oct 2020 at 18:36,
          Nicolas Cadieux <<a
            href="mailto:njacadieux.gitlab@gmail.com"
            moz-do-not-send="true">njacadieux.gitlab@gmail.com</a>>
          wrote:<br>
        </div>
        <blockquote class="gmail_quote" style="margin:0px 0px 0px
          0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">New
          info!<br>
          <br>
          GetProjection works with gdal python 2.2.2 and 2.3.3 but not
          with with <br>
          gdal 3.0.2.  Is this a change in the gdal library?<br>
          <br>
          Nicolas<br>
          <br>
          On 2020-10-26 1:53 p.m., Nicolas Cadieux wrote:<br>
          > Hi,<br>
          ><br>
          > In case mail question was not clear, I have posted the
          questions about <br>
          > this on stackexchange.  I posted a full code that you can
          run exposing <br>
          > my problem and question.<br>
          ><br>
          > Thanks<br>
          ><br>
          > Nicolas<br>
          ><br>
          > <a
href="https://gis.stackexchange.com/questions/377567/cannot-get-projection-from-raster-dataset-using-getprojection"
            rel="noreferrer" target="_blank" moz-do-not-send="true">https://gis.stackexchange.com/questions/377567/cannot-get-projection-from-raster-dataset-using-getprojection</a>
          <br>
          ><br>
          ><br>
          > On 2020-10-25 9:31 p.m., Nicolas Cadieux wrote:<br>
          >> Hi,<br>
          >><br>
          >> With the following code, I get an empty string
          indicating the <br>
          >> projection is not valid.<br>
          >><br>
          >> from osgeo import gdal, osr<br>
          >> raster_ds = gdal.Open(r"C:\temp\180922_WTE3.tif")<br>
          >> target_ds = <br>
          >>
driver.Create(r"c:\temp\output.tif",xsize=raster_ds_ncol,ysize=raster_ds_nrow,bands
          <br>
          >> = 1,eType = gdal.GDT_Float32)<br>
          >> (array is a numpty array.)<br>
          >> raster_srs = osr.SpatialReference()<br>
          >>
          raster_srs.ImportFromWkt(raster_ds.GetProjectionRef())<br>
          >> target_ds.SetProjection(raster_srs.ExportToWkt())<br>
          >> target_ds.GetRasterBand(1).WriteArray(array)<br>
          >> raster_ds = None #close dataset<br>
          >> target_ds=None<br>
          >><br>
          >> Below is the result of gdal info from qgis. File
          appears to have a <br>
          >> valid projection and is properly georeferenced in
          QGIS using other <br>
          >> data sets.   Is this projection wrong or am I missing
          something in my <br>
          >> python code?  The goal is to extract the projection
          from raster_ds <br>
          >> set in order to apply to target_ds.  I can create the
          file, apply a <br>
          >> geotransform but can't get the projection.  Can you
          help?<br>
          >><br>
          >> Thanks/merci!<br>
          >><br>
          >> Nicolas<br>
          >><br>
          >> QGIS version: 3.14.16-Pi<br>
          >> QGIS code revision: df27394552<br>
          >> Qt version: 5.11.2<br>
          >> GDAL version: 3.0.4<br>
          >> GEOS version: 3.8.1-CAPI-1.13.3<br>
          >> PROJ version: Rel. 6.3.2, May 1st, 2020<br>
          >> Processing algorithm…<br>
          >> Algorithm 'Raster information' starting…<br>
          >> Input parameters:<br>
          >> { 'EXTRA' : '', 'INPUT' : 'C:/temp/180922_WTE3.tif',
          'MIN_MAX' : <br>
          >> False, 'NOGCP' : False, 'NO_METADATA' : False,
          'OUTPUT' : <br>
          >> 'TEMPORARY_OUTPUT', 'STATS' : False }<br>
          >><br>
          >> GDAL command:<br>
          >> gdalinfo C:/temp/180922_WTE3.tif<br>
          >> GDAL command output:<br>
          >> Warning 1: TIFFReadDirectory:Sum of Photometric
          type-related color <br>
          >> channels and ExtraSamples doesn't match
          SamplesPerPixel. Defining <br>
          >> non-color channels as ExtraSamples.<br>
          >><br>
          >> Driver: GTiff/GeoTIFF<br>
          >> Files: C:/temp/180922_WTE3.tif<br>
          >> Size is 1194, 2783<br>
          >> Coordinate System is:<br>
          >> PROJCRS["WGS 84 / UTM zone 18N",<br>
          >> BASEGEOGCRS["WGS 84",<br>
          >> DATUM["World Geodetic System 1984",<br>
          >> ELLIPSOID["WGS 84",6378137,298.257223563,<br>
          >> LENGTHUNIT["metre",1]]],<br>
          >> PRIMEM["Greenwich",0,<br>
          >> ANGLEUNIT["degree",0.0174532925199433]],<br>
          >> ID["EPSG",4326]],<br>
          >> CONVERSION["UTM zone 18N",<br>
          >> METHOD["Transverse Mercator",<br>
          >> ID["EPSG",9807]],<br>
          >> PARAMETER["Latitude of natural origin",0,<br>
          >> ANGLEUNIT["degree",0.0174532925199433],<br>
          >> ID["EPSG",8801]],<br>
          >> PARAMETER["Longitude of natural origin",-75,<br>
          >> ANGLEUNIT["degree",0.0174532925199433],<br>
          >> ID["EPSG",8802]],<br>
          >> PARAMETER["Scale factor at natural origin",0.9996,<br>
          >> SCALEUNIT["unity",1],<br>
          >> ID["EPSG",8805]],<br>
          >> PARAMETER["False easting",500000,<br>
          >> LENGTHUNIT["metre",1],<br>
          >> ID["EPSG",8806]],<br>
          >> PARAMETER["False northing",0,<br>
          >> LENGTHUNIT["metre",1],<br>
          >> ID["EPSG",8807]]],<br>
          >> CS[Cartesian,2],<br>
          >> AXIS["(E)",east,<br>
          >> ORDER[1],<br>
          >> LENGTHUNIT["metre",1]],<br>
          >> AXIS["(N)",north,<br>
          >> ORDER[2],<br>
          >> LENGTHUNIT["metre",1]],<br>
          >> USAGE[<br>
          >> SCOPE["unknown"],<br>
          >> AREA["World - N hemisphere - 78°W to 72°W - by
          country"],<br>
          >> BBOX[0,-78,84,-72]],<br>
          >> ID["EPSG",32618]]<br>
          >> Data axis to CRS axis mapping: 1,2<br>
          >> Origin =
          (459417.200000000011642,5028584.700000000186265)<br>
          >> Pixel Size = (0.050000000000000,-0.050000000000000)<br>
          >> Metadata:<br>
          >> AREA_OR_POINT=Area<br>
          >> TIFFTAG_XRESOLUTION=1<br>
          >> TIFFTAG_YRESOLUTION=1<br>
          >> Image Structure Metadata:<br>
          >> INTERLEAVE=BAND<br>
          >> Corner Coordinates:<br>
          >> Upper Left ( 459417.200, 5028584.700) ( 75d31'
          7.03"W, 45d24'34.58"N)<br>
          >> Lower Left ( 459417.200, 5028445.550) ( 75d31'
          6.99"W, 45d24'30.07"N)<br>
          >> Upper Right ( 459476.900, 5028584.700) ( 75d31'
          4.28"W, 45d24'34.59"N)<br>
          >> Lower Right ( 459476.900, 5028445.550) ( 75d31'
          4.24"W, 45d24'30.08"N)<br>
          >> Center ( 459447.050, 5028515.125) ( 75d31' 5.63"W,
          45d24'32.33"N)<br>
          >> Band 1 Block=1194x1 Type=UInt16, ColorInterp=Red<br>
          >> ...<br>
          >> Band 288 Block=1194x1 Type=UInt16,
          ColorInterp=Undefined<br>
          >><br>
          >><br>
          >><br>
          _______________________________________________<br>
          gdal-dev mailing list<br>
          <a href="mailto:gdal-dev@lists.osgeo.org" target="_blank"
            moz-do-not-send="true">gdal-dev@lists.osgeo.org</a><br>
          <a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev"
            rel="noreferrer" target="_blank" moz-do-not-send="true">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a></blockquote>
      </div>
    </blockquote>
  </body>
</html>