<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi,</p>
<p>Looks like it's an environment thing... Both versions of the code
work with 3.1.4. <br>
</p>
<p>Thanks for your help!</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>