Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

So i am trying to use the difference function and i pretty much create the new file i want but the Geometry is empty. No error is produced. any help would great! Thanks.

#!/usr/bin/env python

import os
import sys
import functions
baseFile = 'intlroads.tab'
eraseFile = 'AOI_Motorway.tab'
#performs erase Geoprocessing
#def erase(inputFile, eraseFile): 

try:
  from osgeo import ogr
except ImportError:
   import ogr

try:
   from osgeo import gdal
except ImportError:
  import gdal

driver = ogr.GetDriverByName('Mapinfo File')

feat1 = driver.Open(baseFile, 0)
if feat1 is None:
  print 'Could not open feat1.tab'
  sys.exit(1)

feat2 = driver.Open(eraseFile, 0)
if feat2 is None:
  print 'Could not open feat2.tab'
  sys.exit(1)  


#get layers and loop through them to get geometrys
feat1Layer = feat1.GetLayer()
feat2Layer = feat2.GetLayer()
for feature1 in feat1Layer:
  geomfeat1 = feature1.GetGeometryRef()
for feature2 in feat2Layer:  
  geomfeat2 = feature2.GetGeometryRef()

#Run difference method (matches ESRI's Erase)
GeomDiff = geomfeat1.Difference(geomfeat2)



fn = 'Roads.tab'
if os.path.exists(fn):
  driver.DeleteDataSource(fn)
outDS = driver.CreateDataSource(fn)
if outDS is None:
  print 'Could not create file'
  sys.exit(1)
outLayer = outDS.CreateLayer('roads', geom_type=ogr.wkbLineString)


# copy the fields from the input layer to the output layer
functions.copyFields(feat1Layer, outLayer)



# loop through the input features
inFeature = feat1Layer.GetNextFeature()
while inFeature:
    outLayer = outDS.CreateLayer('roads', geom_type=ogr.wkbLineString)
    #create a new feature
    outFeature = ogr.Feature(featureDefn)

    # set the geometry
    outFeature.SetGeometry(GeomDiff)

    # copy the attributes
    copyAttributes(inFeature, outFeature)

    # add the feature to the output layer
    outLayer.CreateFeature(outFeature)

    # destroy the output feature
    outFeature.Destroy()

    # destroy the input feature and get a new one
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()
share|improve this question

1 Answer 1

up vote 3 down vote accepted

Early in your program, you've done this:

for feature1 in feat1Layer:
    geomfeat1 = feature1.GetGeometryRef()

After which the feat1Layer "cursor" is at the end of the file. Later, when you call

inFeature = feat1Layer.GetNextFeature()

GetNextFeature() will return None. And if it wasn't None, the program would fail at the end of the loop with a NameError due to 'inLayer' being unassigned. If you want to loop over feat1Layer1 twice, you'll either need to close and reopen it, or ResetReading() on it.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.