A bit of a long winded explanation here so please bear with me. I have received grid square data in ETRS89(3035) which I need to cross reference with out own grid squares which are in TM75(29903). However I am having issues with the transformation. I have created some images in QGIS to demonstrate.

The following shows a 10km fishnet created in 3035 bounding Ireland. The pink squares are the 3035 data which I received.

enter image description here

However if I simply re-project my 29903 squares to 3035 you end up with the following issue. You can see that the green squares are tilted. The layer is re-projected but since the squares were drawn using the parallel lines on 29903 they appear tilted in 3035. I realised the issue was likely the same as explained here Projection Transformation from ETRS89 LAEA to WGS1984?. The re-projection is done but the squares need to be essentially redrawn in 3035.

enter image description here

I surmised that the best way to approach this would be to get the centroid of each of my 29903 squares, re-project that point to 3035 and then rebuild my 10km square buffer. I am using shapely to rebuild the buffer. This works..sort of. But as the following image shows the squares are essentially rotated individually and so are slightly out of position or overlapping.

enter image description here

At this point, the way forward which I am considering is iterating through my re-projected centroid points and finding the 3035 square from the 3035 fishnet which contains said centroid and then setting my updated geometry to that square. Does that seem like a valid approach?

The goal is to reproject my 29903 squares to 3035 and have them be on the same grid square scheme as the 3035 data which I need to cross reference. The same applies vice versa in fact. Included below is the script which I am testing to do the reprojection. Note I have not implemented my fishnet approach yet.

import ogr, osr
import os
from shapely.geometry import Point

driver = ogr.GetDriverByName('ESRI Shapefile')

# input SpatialReference
inSpatialRef = osr.SpatialReference()

# output SpatialReference
outSpatialRef = osr.SpatialReference()

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

# get the input layer
inDataSet = driver.Open(r'I:\Code\Python\input\test_shape\test_shape.shp')
inLayer = inDataSet.GetLayer()

# create the output layer
outputShapefile = r'I:\Code\Python\output\test_shape\test_shape.shp'
if os.path.exists(outputShapefile):
outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer("basemap", outSpatialRef, geom_type=ogr.wkbMultiPolygon)

# add fields
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)

# 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()
    centroid = geom.Centroid()
    geom = Point(centroid.GetX(), centroid.GetY())
    geom = geom.buffer(5000, cap_style=3)
    geom = ogr.CreateGeometryFromWkt(geom.wkt)

    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    # dereference the features and get the next input feature
    outFeature = None
    inFeature = inLayer.GetNextFeature()

# Save and close the shapefiles
inDataSet = None
outDataSet = None

Your Answer

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

Browse other questions tagged or ask your own question.