Majelika

Reprojecting a shapefile with GDAL OGR and Python

I needed to reproject shapefiles -- unprojected lat/long WGS84 -- and once again the page Geoprocessing with Python using Open Source GIS was of great help, so the credit goes to this page, which is very worthwhile studying as an introduction to using gdal/ogr and Python.

The code below is an adjustment of homework 2b in the page above for my case.

I tried to add explanations to the original code so that its easy to follow by just reading it -- hope Im right....

# -*- coding: utf-8 -*-  
# Import Modules 
import ogr, osr, os, sys 
 
# set file names 
infile = C:UsersmaxDocumentsIcechartsDataIceChart_2012ice20120608.shp 
outfile = C:UsersmaxDocumentsIcechartsDataIceChart_2012ice20120608_epsg3575.shp 
 
#get path and filename seperately 
(outfilepath, outfilename) = os.path.split(outfile) 
#get file name without extension            
(outfileshortname, extension) = os.path.splitext(outfilename)   #get file name without extension 
 
# Spatial Reference of the input file 
# Access the Spatial Reference and assign the input projection 
inSpatialRef = osr.SpatialReference() 
inSpatialRef.ImportFromEPSG(4326)         # unprojected WGS84 
 
# Spatial Reference of the output file 
# Access the Spatial Reference and assign the output projection 
# UTM 33N which we use for all Spitsbergen does not work since area too far  
# outside UTM 33 
outSpatialRef = osr.SpatialReference() 
 outSpatialRef.ImportFromEPSG(3575) # EPSG:3575: WGS 84 / North Pole LAEA Europe 
  
# create Coordinate Transformation 
coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef) 
 
# Open the input shapefile and get the layer 
driver = ogr.GetDriverByName(ESRI Shapefile) 
indataset = driver.Open(infile, 0) 
if indataset is None: 
    print  Could not open file 
    sys.exit(1) 
inlayer = indataset.GetLayer() 
 
# Create the output shapefile but check first if file exists 
if os.path.exists(outfile): 
    driver.DeleteDataSource(outfile) 
 
outdataset = driver.CreateDataSource(outfile) 
 
if outfile is None: 
    print  Could not create file 
    sys.exit(1) 
outlayer = outdataset.CreateLayer(outfileshortname, geom_type=ogr.wkbPolygon) 
 
# Get the FieldDefn for attributes ICE_TYPE and ID and add to output shapefile 
# (which i know are in my file) 
feature = inlayer.GetFeature(0) 
fieldDefn1 = feature.GetFieldDefnRef(ICE_TYPE) 
fieldDefn2 = feature.GetFieldDefnRef(ID) 
 
outlayer.CreateField(fieldDefn1) 
outlayer.CreateField(fieldDefn2) 
 
# get the FeatureDefn for the output shapefile 
featureDefn = outlayer.GetLayerDefn() 
 
#Loop through input features and write to output file 
infeature = inlayer.GetNextFeature() 
while infeature: 
  
    #get the input geometry 
    geometry = infeature.GetGeometryRef() 
  
    #reproject the geometry, each one has to be projected seperately 
    geometry.Transform(coordTransform) 
  
    #create a new output feature 
    outfeature = ogr.Feature(featureDefn) 
  
    #set the geometry and attribute 
    outfeature.SetGeometry(geometry) 
    outfeature.SetField(ICE_TYPE, infeature.GetField(ICE_TYPE)) 
    outfeature.SetField(ID, infeature.GetField(ID)) 
  
    #add the feature to the output shapefile 
    outlayer.CreateFeature(outfeature) 
  
    #destroy the features and get the next input features 
    outfeature.Destroy 
    infeature.Destroy 
    infeature = inlayer.GetNextFeature() 
  
#close the shapefiles 
indataset.Destroy() 
outdataset.Destroy() 
 
#create the prj projection file 
outSpatialRef.MorphToESRI() 
file = open(outfilepath + \+ outfileshortname + .prj, w) 
file.write(outSpatialRef.ExportToWkt()) 
file.close() 
 
#end





The unprojected file:


The projected shape file:

Posted by Unknown at 3:45 PM

Email ThisBlogThis!Share to TwitterShare to FacebookShare to Pinterest

Labels: a, and, gdal, ogr, python, reprojecting, shapefile, with