[gdal-dev] GetProjection returns empty string, python 3.7
Paul Harwood
runette at gmail.com
Mon Oct 26 16:18:35 PDT 2020
"and the Anaconda Navigator does not give me higher options"
just run the command "conda install -c conda-forge gdal".
It has to come from conda-forge and not conda-main. You can add conda-forge
to the navigator - but it is usually just quicker to use the cli.
Sorry - but I don't really have the time to debug test files.
On Mon, 26 Oct 2020 at 23:00, Nicolas Cadieux <njacadieux.gitlab at gmail.com>
wrote:
> Hi,
>
> 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?
>
> Thanks, I appreciate the 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.gitlab at 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.gitlab at 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.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
>>>> 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/84a54dac/attachment-0001.html>
More information about the gdal-dev
mailing list