Re: How can I "clip" shapefiles to a box-shaped subset without using the ESRI ArcMap tool?

Bryan Keith <[email protected]> Mon, 28 Aug 2006 08:35:45 -0600

Newsgroups	gmane.comp.gis.freegis
Message-ID	<[email protected]>

Murray,

I started writing a Python script using ogr to do this about 10 months 
ago.  As I recall, it worked for what I needed at the time, but it's 
still pretty buggy.  I have only one version in CVS and another version 
that looks hacked up.  I attached the CVS version.  Be warned that I 
haven't looked at this in quite some time and recall that it needs some 
work.  Maybe it will be a useful place for you to start.

Bryan

Bryan Keith
GIS Specialist
Geomega, Inc.
Boulder, CO, USA

Murray Brown wrote:
> Our marine data training program (www.oceanteacher.org) tries to
> avoid necessary use of commercial software, and so far we've been
> able to stay nearly ESRI-free.  But the one Arc tool that we still
> find necessary is the clip function that cuts a shapefile down to the
> borders defined by a second shapefile.  In our case, the second
> shapefile is almost always a simple rectangular box (area of
> interest).
> 
> Does anybody know a free, stand-alone program (hopefully a binary
> executable) that can do this also?
> 
> Thanks in advance for any help.
> 
> Murray Brown UNESCO-IOC Marine Data Mgt. Training Program
> 
> 
> ------------------------------------------------------------------------
> 
> 
> _______________________________________________ Freegis-list mailing
> list [email protected] 
> https://intevation.de/mailman/listinfo/freegis-list

_______________________________________________
Freegis-list mailing list
[email protected]
https://intevation.de/mailman/listinfo/freegis-list

clip.py (text/plain, 2.7 KB)

#This routine clips a shapefile.
#Original coding Oct 18, 2005 by Bryan
#currently does not work on win32 because GEOS isn't supported in FWTools on win32
#(true as of version 1.0.0a4)
import ogr
import os
import sys
import string
import platform
def Usage():
	print "Usage: clip inshapefile wktPYfile outshapefile"
	return
	
def AddFeature(ogrSLayerNew, newGeom, myFeature):
	newFeature = ogr.Feature(ogrSLayerNew.GetLayerDefn())
	newFeature.SetGeometry(newGeom)
	for k in range(myFeature.GetFieldCount()):
		newFeature.SetField(k,myFeature.GetField(k))
		
	ogrSLayerNew.CreateFeature(newFeature)
				
if len(sys.argv) != 4:
	Usage()
	sys.exit()

#check that we're on Linux
if platform.system() != "Linux":
	print "This tool is currently only available on Linux because GEOS isn't " + \
		"supported on "
	print "win32."
	sys.exit()
	
#remove .shp from outshapefile
strNewFN = sys.argv[3]
if (len(strNewFN) > 4):
	if string.lower(strNewFN[-4:]) == ".shp":
		strNewFN = strNewFN[:-4]
		
		
#remove .shp from inshapefile
strFN = sys.argv[1]
if (len(strFN) > 4):
	if string.lower(strFN[-4:]) == ".shp":
		strFN = strFN[:-4]		
		
		
#check if input file exists
if not os.path.isfile(strFN + ".shp"):
	print strFN + ".shp does not exist"
	sys.exit()
	
#check if output file exists
if os.path.isfile(strNewFN + ".shp"):
	print strNewFN + ".shp exists"
	sys.exit()

#check if wktPYfile exists
if not os.path.isfile(sys.argv[2]):
	print sys.argv[2] + " does not exist"
	sys.exit()
	
ogrSDS = ogr.Open(os.path.abspath(""),1)
ogrSLayer = ogrSDS.GetLayerByName(strFN)
myFeatureDefn = ogrSLayer.GetLayerDefn()
ogrSLayerNew = ogrSDS.CreateLayer(strNewFN, ogrSLayer.GetSpatialRef(), \
	myFeatureDefn.GetGeomType())
#add fields
for i in range(myFeatureDefn.GetFieldCount()):
	myFieldDefn = myFeatureDefn.GetFieldDefn(i)
	ogrSLayerNew.CreateField(myFieldDefn)
	
#make clipping geometry
wkt = open(sys.argv[2]).read()
clipGeom = ogr.CreateGeometryFromWkt(wkt)
#loop through the input features
for i in range(ogrSLayer.GetFeatureCount()):
	myFeature = ogrSLayer.GetFeature(i)
	oldGeom = myFeature.GetGeometryRef()
	newGeomBag = clipGeom.Intersection(oldGeom)
	#newGeomBag = oldGeom.Intersection(clipGeom)
	if newGeomBag is not None:
		#how do I loop through each geometry?
		#print str(newGeomBag.GetGeometryCount())
		if newGeomBag.GetGeometryCount() > 0:
			for j in range(newGeomBag.GetGeometryCount()):
				newGeom = newGeomBag.GetGeometryRef(j)
				#add this feature
				AddFeature(ogrSLayerNew, newGeom, myFeature)
								
		else:
			#add this feature
			if newGeomBag.GetGeometryType() == oldGeom.GetGeometryType():
				AddFeature(ogrSLayerNew, newGeomBag, myFeature)
				
				
		
			
ogrSLayerNew.SyncToDisk()
ogrSDS.Destroy()
print strFN + ".shp clipped with " + sys.argv[2] + " to create " + strNewFN + ".shp"