0

I'm trying to get the location coordinates written to a shp file after reprojecting it. I have shp files that are in WGS 84 (EPSG:4326) and I reproject them into UTM 13N (EPSG:27700). I then want to add two columns and populate them with the UTM coordinates, the new projection. I know that this is easy to do with QGIS in the attribute table but I have several files and would like to find a way to batch the process and handle several files at once so I am hoping there is an OGR method.

1

You can try the following script and loop through every shapefile you need to reproject. Alternatively, ogr will treat a folder full of shapefiles as a data set (ds) and each shapefile as a layer (lyr), so you could simply declare ds as the folder where all your shapefiles are and loop through the layers (shapefiles).

from osgeo import ogr, osr

# Create the UTM Spatial Reference
sr = osr.SpatialReference()
sr.ImportFromEPSG(27700)

# Read existing shapefile
fn = "shapefile"
ds = ogr.Open(fn)
lyr = ds.GetLayer()

# Create new shapefile
new_fn = "new_shapefile"
driver = ogr.GetDriverByName('ESRI Shapefile')
new_ds = driver.CreateDataSource(new_fn)
new_lyr = new_ds.CreateLayer('UTM', sr, ogr.wkbPoint)

# Copy the original attrbiute fields
new_lyr.CreteFields(lyr.schema)

# Create the two UTM attribute fields
coord_fld = ogr.FieldDefn('UTM_x', ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
new_lyr.CreateField(coord_fld)
coord_fld.SetName('UTM_y')
new_lyr.CreateField(coord_fld)

# Create dummy feature (to re-use) for the new shapefile
new_feat = ogr.Feature(new_lyr.GetLayerDefn())

# Loop through the features in the original shapefile
for feat in lyr:

    # Get input geometry
    geom = feat.geometry().Clone()

    # Transform and set new geometry
    geom.TransformTo(sr)
    new_feat.SetGeometry(geom)

    # Populate original attribute fields
    for i in range(feat.GetFieldCount()):
        value = feat.GetField(i)
        new_feat.SetField(i, value)

    # Populate UTM attribute fields
    new_feat.SetField('UTM_x', geom.GetX())
    new_feat.SetField('UTM_y', geom.GetY())

    # Create actual new feature
    new_lyr.CreateFeature(new_feat)

del ds, new_ds
  • @Marcello Villa Thanks for the direction. However, I can't get this to work. I made fn the path to a .shp file – Casivio Feb 9 '18 at 19:15
  • @Casivio Could you elaborate why it is not working? – Marcelo Villa Feb 9 '18 at 19:21
  • @Marcello Villa I was and accidentally pressed enter so it posted. I set fn = 'path to a .shp file' and new_fn = 'path with a new name'. This produces the error: new_lyr.CreateField(coord_fld) AttributeError: 'NoneType' object has no attribute 'CreateField'. I also tried setting ds = ogr.Open('path to folder with .shp files) and it produced this error: lyr = ds.GetLayer() AttributeError: 'NoneType' object has no attribute 'GetLayer' I tried ds = 'path to folder' which did not work because it was a string. I am using Pythong 2.7 in case that makes a difference – Casivio Feb 9 '18 at 19:26
  • @Casivio It looks like ogr is not being able to open the shapefile nor create the new one. Are yopu sure the paths are correct? Could you paste them here? – Marcelo Villa Feb 9 '18 at 22:44
  • @Marcello Villa I do not think that is the case. I'm using fn = path+os.sep+'Data'+os.sep+'shp Location Files'+os.sep+'ESRI Shapefile'+os.sep+'WGS_84'+os.sep+'UTM_13N_Americium-241.shp' which prints out as 'C:\\Projects\\770- LANL\\SLVs\\Data\\shp Location Files\\ESRI Shapefile\\WGS_84\\UTM_13N_Americium-241.shp' when called. – Casivio Feb 9 '18 at 23:20

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.