1

I have 4 shapefiles in the same folder names like that (1996.shp, 2000.shp, 2006.shp, 2014.shp) and one raster.

I need to reproject all shapfiles in the same projection (4326) of raster, but I do it for one shapefile that it script :

from osgeo import ogr, osr, gdal

    #tif with projections I want
    tif = gdal.Open("D:/projetp/NDVI.tif")

    #shapefile with the from projection
    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataSource =   driver.Open("C:/projetp/shape/1996.shp", 1)
    layer = dataSource.GetLayer()

    #set spatial reference and transformation
    sourceprj = layer.GetSpatialRef()
    targetprj = osr.SpatialReference(wkt = tif.GetProjection())
    transform = osr.CoordinateTransformation(sourceprj, targetprj)

    to_fill = ogr.GetDriverByName("Esri Shapefile")
    ds = to_fill.CreateDataSource("C:/projetp/shape/1996poly.shp")
    outlayer = ds.CreateLayer('', targetprj, ogr.wkbPolygon)
    outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))

    #apply transformation

    i = 0
    for feature in layer:
        transformed = feature.GetGeometryRef()
        transformed.Transform(transform)

        geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb())
        defn = outlayer.GetLayerDefn()
        feat = ogr.Feature(defn)
        feat.SetField('id', i)
        feat.SetGeometry(geom)
        outlayer.CreateFeature(feat)
        i += 1
        feat = None

    ds = None

How can I execute the script for all shapefiles in the folder?

I'm new in Python 3.

Improve this question
1

You could wrap your code into one function (e.g. reprojectShape) and then call it from another function which reads files within folder where shapefiles are located. Here is the code:

from osgeo import ogr, osr, gdal
import os

def reprojectShape(fullPath, filename):
  #tif with projections I want
  tif = gdal.Open("D:/projetp/NDVI.tif")

  #shapefile with the from projection
  driver = ogr.GetDriverByName("ESRI Shapefile")
  dataSource =   driver.Open(fullPath + filename, 1)
  layer = dataSource.GetLayer()

  #set spatial reference and transformation
  sourceprj = layer.GetSpatialRef()
  targetprj = osr.SpatialReference(wkt = tif.GetProjection())
  transform = osr.CoordinateTransformation(sourceprj, targetprj)

  to_fill = ogr.GetDriverByName("Esri Shapefile")
  ds = to_fill.CreateDataSource(fullPath + 'reprojected/' + filename)
  outlayer = ds.CreateLayer('', targetprj, ogr.wkbPolygon)
  outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))

  #apply transformation

  i = 0
  for feature in layer:
    transformed = feature.GetGeometryRef()
    transformed.Transform(transform)

    geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb())
    defn = outlayer.GetLayerDefn()
    feat = ogr.Feature(defn)
    feat.SetField('id', i)
    feat.SetGeometry(geom)
    outlayer.CreateFeature(feat)
    i += 1
    feat = None

  ds = None

def reprojectAll():
  folderShape = "C:/projetp/shape/"
  for filename in os.listdir(folderShape):
    if filename.endswith('.shp'):
      reprojectShape(folderShape, filename)
Improve this answer
2
  • Vranic thank you , but in the fullpath and filname what i need to do ? – Majda El yacouby Jan 18 '19 at 14:10
  • def reprojectShape(fullPath, filename): #tif with projections I want tif = gdal.Open("D:/projetp/NDVI.tif") #shapefile with the from projection driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open(fullPath + filename, 1) layer = dataSource.GetLayer() need to defined path like that path=('C:/projetp/shape/1996.shp") and filename=("1996.shp") ? – Majda El yacouby Jan 18 '19 at 14:43
1

Note that there are tools to directly reproject your shapefiles (not feature by feature), e.g. ogr2ogr

import subprocess, glob
ogr2ogrpath = "C:\\Users\\maj\\Anaconda3\\Library\bin\\ogr2ogr.exe" #local path provided in the comments
for shp in glob.glob("D:/projetp/shape/*shp"):#list all shapefiles in a folder
    subprocess.call([ogr2ogrpath, "-t_srs", "EPSG:4326", shp.replace("projetp/shape","reprojected"), shp]) #remark, you might need to give full path to your ogr2ogr.exe (if ogr2ogr.exe isnot in your PATH list)
Improve this answer
10
  • thank you i can use you script to do the same work to my script or i need add your script ? – Majda El yacouby Jan 18 '19 at 9:36
  • 1
    my script will do the same job as your script, except that I do not select the coordiante system dynamically based on a given image (EPSG:4326 is hardcoded). You can use what you did in your script to get your targetprj and replace my "EPSG:4326" with targetprj if you need a dynamic projection – radouxju Jan 18 '19 at 10:42
  • ` #remark, you might need to give full path to your ogr2ogr.exe (if ogr2ogr.exe isnot in your PATH list)` i don't understand what i need to do in path – Majda El yacouby Jan 18 '19 at 13:59
  • I am not using the API but I use subprocess to run a command line. If you have the API, you must have the ogr2ogr binary file installed somewhere on you computer. If you installed ogr2OGR using OSGEO4W, for instance, run this script from within the OSGEO4W command prompt will work. Otherwise, you must find where it is located and write its full path (e.g. C:/osgeo/bin/ogr2ogr.exe) – radouxju Jan 18 '19 at 14:50
  • i'm using spyder anaconda i don't know how do that – Majda El yacouby Jan 18 '19 at 17:37

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.