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.git...@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.git...@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.git...@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@lists.osgeo.org >>> https://lists.osgeo.org/mailman/listinfo/gdal-dev >> >>
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev