<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
  </head>
  <body>
    <p>Hi,</p>
    <p>I am currently on 3.0.2.  and the Anaconda Navigator does not
      give me higher options.  I will create a new clean environments
      and try conda install.  I will try changing the code I think my
      version of the code probably corresponds to <3 versions of
      gdal.  It's a cut and paste from the gdal tutorial. Can I send you
      my test file to test on your code?  <br>
    </p>
    <p>Thanks, I appreciate the help.<br>
    </p>
    <p>Nicolas<br>
    </p>
    <div class="moz-cite-prefix">On 2020-10-26 6:31 p.m., Paul Harwood
      wrote:<br>
    </div>
    <blockquote type="cite"
cite="mid:CAE8nN5P4o5CNh6vDrksN=fCUBTR0q3LMjYWEGP9v11A_dS-08g@mail.gmail.com">
      <meta http-equiv="content-type" content="text/html; charset=UTF-8">
      <div dir="ltr">Me too.
        <div><br>
        </div>
        <div>But which conda package for gdal are you using - there are
          a lot of ones out there. </div>
        <div><br>
        </div>
        <div>Are you using conda install -c conda-forge gdal ? </div>
        <div><br>
        </div>
        <div>It should be version 3.1.4 at the moment ?</div>
        <div><br>
        </div>
        <div>Also - I don't know if this is part of it but ...</div>
        <div><br>
        </div>
        <div>"srs = osr.SpatialReference()</div>
        srs.SetUTM(11, 1)<br>
        srs.SetWellKnownGeogCS("NAD27")<br>
        dst_ds.SetProjection(srs.ExportToWkt())
        <div>"</div>
        <div><br>
        </div>
        <div>why not just</div>
        <div><br>
        </div>
        <div>
          <div>"srs = osr.SpatialReference()</div>
          srs.SetUTM(11, 1)<br>
          srs.SetWellKnownGeogCS("NAD27")<br>
          dst_ds.SetSpatialRef(srs)
          <div>"</div>
          <div>??</div>
        </div>
      </div>
      <br>
      <div class="gmail_quote">
        <div dir="ltr" class="gmail_attr">On Mon, 26 Oct 2020 at 22:11,
          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">
          <div>
            <p>Hi,</p>
            <p>Comes back as None for me... :(  I'am working with
              Anaconda.  Could it be my environment?</p>
            <p>Nicolas<br>
            </p>
            <div>On 2020-10-26 5:04 p.m., Paul Harwood wrote:<br>
            </div>
            <blockquote type="cite">
              <div dir="ltr">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.
                <div><br>
                </div>
                <div>This code is working for me at this very minute and
                  giving me the name of the CRS for a gtiff</div>
                <div><br>
                </div>
                <div>dataset = gdal.Open(str(filepath))<br>
                          if dataset:<br>
                              proj = "None"<br>
                              crs = dataset.GetSpatialRef()<br>
                              if crs:<br>
                                  proj = crs.GetName()<br>
                              return {"type": "gdal", "driver":
                  dataset.GetDriver().ShortName, "proj": proj}<br>
                </div>
              </div>
              <br>
              <div class="gmail_quote">
                <div dir="ltr" class="gmail_attr">On Mon, 26 Oct 2020 at
                  20:08, Nicolas Cadieux <<a
                    href="mailto:njacadieux.gitlab@gmail.com"
                    target="_blank" 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">
                  <div>
                    <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
                        href="https://gdal.org/tutorials/raster_api_tut.html"
                        target="_blank" moz-do-not-send="true">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
href="https://stackoverflow.com/questions/37648439/simplest-way-to-save-array-into-raster-file-in-python"
                        target="_blank" moz-do-not-send="true">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
                        href="https://gdal.org/tutorials/raster_api_tut.html"
                        target="_blank" moz-do-not-send="true">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
href="https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/gdalinfo.py"
                        target="_blank" moz-do-not-send="true">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>On 2020-10-26 3:44 p.m., Paul Harwood wrote:<br>
                    </div>
                    <blockquote type="cite">
                      <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"
                            target="_blank" 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>
                  </div>
                </blockquote>
              </div>
            </blockquote>
          </div>
        </blockquote>
      </div>
    </blockquote>
  </body>
</html>