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.
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.
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.
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() inSpatialRef.ImportFromEPSG(29903) # output SpatialReference outSpatialRef = osr.SpatialReference() outSpatialRef.ImportFromEPSG(3035) # 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): driver.DeleteDataSource(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) 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.Centroid().GetX()) centroid = geom.Centroid() centroid.Transform(coordTrans) geom = Point(centroid.GetX(), centroid.GetY()) geom = geom.buffer(5000, cap_style=3) geom = ogr.CreateGeometryFromWkt(geom.wkt) geom.AssignSpatialReference(outSpatialRef) print(geom.Centroid().GetX()) # 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) # dereference the features and get the next input feature outFeature = None inFeature = inLayer.GetNextFeature() # Save and close the shapefiles inDataSet = None outDataSet = None