1

A have input .tif files and geojson with different EPGS. I'm trying to use osr.CoordinateTransformation(source, target) method which takes 2 arguments of EPGS like on this cookbook. So how can I get this values from input files? Method AutoIdentifeEPSG () doesn't work.

Here is my code:

Register Imagine driver and open file. Getting projection of input file.

driver = gdal.GetDriverByName('GTiff')
dataset = gdal.Open('{}/{}.tif'.format(image_dir, file_name))
projection = dataset.GetProjection()

And I created geometry to set it directly to the feature

geometry = ogr.CreateGeometryFromJson(""" some json string """)
lyr = memds.CreateLayer('', geom_type=ogr.wkbMultiPolygon)
feat = ogr.Feature(lyr.GetLayerDefn())
feat.SetGeometryDirectly( ogr.Geometry(wkt = geometry.ExportToWkt()) )
lyr.CreateFeature(feat)

Then I need to transform geomerty projection to the projection of the input .tif file.

geometry.Transform(osr.CoordinateTransformation(…), …)

but when I tried to use .Transform () method on my geometry I got error like: `NoneType object have not attribute or method Transform.

  • In my answer you also have how to use 'AutoIdentifyEPSG' method for getting EPSG code. – xunilk Jul 6 '16 at 18:24
3

Next code is similar to the cookbook code for testing the OGR 'Transform' method (from EPSG:32612 to EPSG:4326). It also includes 'AutoIdentifyEPSG' and 'GetAuthorityCode' methods for determining EPSG code for projection in my GeoJSON layer (EPSG:32612).

from osgeo import ogr, osr
import os

driver = ogr.GetDriverByName('GeoJSON')

# get the input layer
inDataSet = driver.Open("/home/zeito/pyqgis_data/Points.geojson")
inLayer = inDataSet.GetLayer()

inSpatialRef = inLayer.GetSpatialRef()

# loading projection
sr = osr.SpatialReference(str(inSpatialRef))

# detecting EPSG/SRID
res = sr.AutoIdentifyEPSG()

srid = sr.GetAuthorityCode(None)

print srid

# input SpatialReference
inSpatialRef.ImportFromEPSG(int(srid))

# output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# create the output layer
outputShapefile = "/home/zeito/pyqgis_data/Points_4326.geojson"

if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)

outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer("Points_4326", geom_type=ogr.wkbPoint)

# add fields
inLayerDefn = inLayer.GetLayerDefn()

for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
inFeature = inLayer.GetNextFeature()

while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    print geom
    # reproject the geometry
    geom.Transform(coordTrans)
    print geom
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # destroy the features and get the next input feature
    outFeature.Destroy()
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()

# close the shapefiles
inDataSet.Destroy()
outDataSet.Destroy()

After running the code at the Python Console of QGIS (see next image), it can be observed the ESPG code of the layer (it was adequately imported from GeoJSON layer) and the geometry for each feature (WKT format) before and after 'Transform' method application.

enter image description here

It can be also observed that the GeoJSON output layer, loaded at the Map Canvas of QGIS, was adequately reprojected from UTM12N/WGS84 projection.

enter image description here

3

To transform from lat/lon (4326) to any coordinate system I do following:

# Getting spatial reference of input raster
srs = osr.SpatialReference()
srs.ImportFromWkt(projection)

# WGS84 projection reference
OSR_WGS84_REF = osr.SpatialReference()
OSR_WGS84_REF.ImportFromEPSG(4326)

# OSR transformation
wgs84_to_image_trasformation = osr.CoordinateTransformation(OSR_WGS84_REF, srs)

Create geometry. Finally:

wkt_geom.Transform(wgs84_to_image_trasformation)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.