arcgis-desktoparcmaperase
I'm wondering if there's a way to intersect two polygon layers, take the intersection, and remove that from one of the polygon layers, making a new layer out of whatever is left over?
Basically, the idea I have is one of "punching holes" in a map of Africa where the holes are points around which I construct some sort of buffer.
Is there a built-in or user created script for this using ArcGIS Desktop 10.0?
I'm unclear on the drawing portion of your question, so I hope I'm not totally off base here, but it seems to me that you could apply a buffer based on the z-value field around each of your points. Take the resulting polygon and use the polygon to raster tool to create your final raster (if you have ArcInfo). If you don't have ArcInfo I believe gdal_rasterize should do the trick.
Note: If your buffer distance around each point isn't directly correlated to your z-value, just add an extra buffer distance field.
ESRI link: http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=Polygon_to_Raster_(conversion)
gdal_rasterize link: http://www.gdal.org/gdal_rasterize.html
I am using vector and raster approach to solve it. The script below has 3 input parameters: a) polygons (must be single part), b)output point feature class to store the points most distant from polygon boundary c) tolerance in map units.
Algorithm:
POLYGON =shape geometry
The idea borrowed if memory serves from very old post by @whuber.
import arcpy, traceback, os, sys
from arcpy import env
env.overwriteOutput = True
infc = arcpy.GetParameterAsText(0)
outfc=arcpy.GetParameterAsText(1)
tolerance=float(arcpy.GetParameterAsText(2))
d=arcpy.Describe(infc)
SR=d.spatialReference
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
theFolder,theFile=os.path.split(outfc)
arcpy.CreateFeatureclass_management(theFolder, theFile, "POINT", infc, "DISABLED", "DISABLED", SR)
arcpy.AddField_management(outfc, "theDist", "DOUBLE")
shapefield,fidName = d.ShapeFieldName,d.OIDFieldName
fileds2add=[]
fields = [f.name for f in arcpy.ListFields(infc)]
for f in fields:
if f not in (shapefield,fidName,"Shape_Length","Shape_Area"):
fileds2add.append(f)
tbl=arcpy.da.TableToNumPyArray(infc,fileds2add)
nF=len(tbl)
arcpy.SetProgressor("step", "", 0, nF,1)
arcpy.AddMessage("Creating points... ")
fileds2add=["SHAPE@","theDist"]+fileds2add
curT = arcpy.da.InsertCursor(outfc,fileds2add)
theM=0
with arcpy.da.SearchCursor(infc, "SHAPE@") as rows:
for row in rows:
feat = row[0]; prt=feat.getPart(0)
feat=arcpy.Polygon(prt,SR)
outLine=feat.boundary(); pSave=feat.trueCentroid
d=outLine.distanceTo(pSave)
if d<=tolerance:break
while (True):
theLine=feat.boundary()
p=feat.labelPoint
d=theLine.distanceTo(p)
try:
buf=theLine.buffer(d)
except:
pSave=feat.labelPoint
break
intR=feat.difference(buf)
n=intR.partCount;aMax=0
for i in xrange (n):
prt=intR.getPart(i)
pgon=arcpy.Polygon(prt,SR)
aCur=pgon.area
if aCur>aMax:
aMax=aCur;feat=pgon
pN=feat.labelPoint
d=arcpy.PointGeometry(p).distanceTo(pN)
if d<=tolerance:
pSave=pN; break
d=outLine.distanceTo(pSave)
theRow=[pSave,d]; theP=list(tbl[theM])
theRow+=theP
curT.insertRow(theRow)
theM+=1
arcpy.SetProgressorPosition()
del tbl
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Script adds field theDist to point table as well as all visible attributes of polygon input feature class. Use this distance to make relevant buffer around points, or just select the points that have theDist<=your threshold.
Note it is very slow (10 polygons like one on the picture per second). Try on selection and use reasonable tolerance.
There is faster raster solution:
UPDATE Nov.2016 The below script is modification of original. It handles donut polygons and works slightly faster, approximately 25 polygons/second:
import arcpy, traceback, os, sys
from arcpy import env
env.overwriteOutput = True
infc = arcpy.GetParameterAsText(0)
outfc=arcpy.GetParameterAsText(1)
tolerance=float(arcpy.GetParameterAsText(2))
d=arcpy.Describe(infc)
SR=d.spatialReference
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
theFolder,theFile=os.path.split(outfc)
arcpy.CreateFeatureclass_management(theFolder, theFile, "POINT", infc, "DISABLED", "DISABLED", SR)
arcpy.AddField_management(outfc, "theDist", "DOUBLE")
shapefield,fidName = d.ShapeFieldName,d.OIDFieldName
fileds2add=[]
fields = [f.name for f in arcpy.ListFields(infc)]
for f in fields:
if f not in (shapefield,fidName,"Shape_Length","Shape_Area"):
fileds2add.append(f)
tbl=arcpy.da.TableToNumPyArray(infc,fileds2add)
nF=len(tbl)
arcpy.SetProgressor("step", "", 0, nF,1)
arcpy.AddMessage("Creating points... ")
fileds2add=["SHAPE@","theDist"]+fileds2add
curT = arcpy.da.InsertCursor(outfc,fileds2add)
with arcpy.da.SearchCursor(infc, "SHAPE@") as rows:
for m,row in enumerate(rows):
feat = row[0]
outLine=feat.boundary()
pSave=feat.labelPoint
d=outLine.distanceTo(pSave)
while (True):
theLine=feat.boundary()
p=feat.labelPoint
dist=theLine.distanceTo(p)
try:
buf=feat.buffer(-dist)
except:
pSave=feat.labelPoint
break
try: pN=buf.labelPoint
except:pN=feat.labelPoint
d=arcpy.PointGeometry(p).distanceTo(pN)
if d<=tolerance:
pSave=pN; break
feat=buf
d=outLine.distanceTo(pSave)
theRow=[pSave,d]; theP=list(tbl[m])
theRow+=theP
curT.insertRow(theRow)
arcpy.SetProgressorPosition()
del tbl
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Best Answer
It sounds like you want Symmetric Difference
In ArcGIS you could use the Symmetrical Difference tool http://resources.arcgis.com/en/help/main/10.1/index.html#//00080000000r000000
QGIS has it too: