[gdal-dev] GetProjection returns empty string, python 3.7
Nicolas Cadieux
njacadieux.gitlab at gmail.com
Mon Oct 26 13:08:39 PDT 2020
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.gitlab at gmail.com <mailto:njacadieux.gitlab at 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 at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org>
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201026/369ef445/attachment-0001.html>
More information about the gdal-dev
mailing list