Latest Contributions by XanderBakker

Hi Erin, Thanks for posting the data. I see what you mean, the TIN will not cover the entire area needed by the line. To solve this, I guess some programming would be needed. I started with something in Python, but the possibilities offered in ArcObjects for working with geometries is still superior (for instance snapping a point to the line and determining the M-position). For how many lines would you like to do this? What do you want to do with the end points (should the slope be extrapolated towards the beginning and end based on the slope determined by the two points located nearest to these from and end points? Or should the last value (closest to from or to point) be used for these points? I can look into this a little more, but it is going to take some time and the projects I'm working on currently don't give me much room for this. Kind regards, Xander ... View more

Hi Erin, Is it possible to share a small part of your data? That makes it easier to see what can be done. Kind regards, Xander ... View more

Hi Mark, Basically the initial idea was a search distance of 200m (actually 200-400m, but I wanted to keep the array small). If you have a pixels size of 10m, a distance of 200m will be 20 pixels. So twenty pixels on each side means 20 * 2 = 40. When using filters and similar calculations, one uses an impair number of pixels to have a central pixel (the one where the value will be assigned). This is why the size is 41 pixels (2 * 20 + 1) and the winhalf is 20 pixels. Kind regards, Xander ... View more

Hi Neil, I looked at your original post: http://forums.arcgis.com/threads/97771-Calculate-percentage-of-NODATA-pixels-in-a-scene?p=347836#post347836 ...and the fact that you can't access the RAT. This is strange. I did something similar in this thread: http://forums.arcgis.com/threads/97645-Export-raster-attribute-table-into-CSV-using-Python-automation-! ... I used the da cursor to access the RAT of a list of rasters and it worked. I tried to reproduce your error and in my case it also happened. Looking at the differences I noticed that in my case I don't supply a raster object to the cursor, but the path to the raster. In that case (apparently) arcpy has the possibility to interpret it as a "table" (at least when it's an integer raster with a RAT) ... So change it to: locRas1 = os.path.join(env.workspace,'FlowD') for r in arcpy.da.SearchCursor(locRas1, ["VALUE", "COUNT"]): Then it should work. I think you should mention this to Esri Support, since it would be a good thing if you can use a raster object in a da cursor... Kind regards, Xander ... View more

Hi Anne, I notice a few more things. I wrote the code based on rasters and points FC in a file geodatabase. You are using grids and a point shapefile. One of the errors is "No valid Observation points present.". This is likely the cause of the wrong whereclause that is being build. The whereclause is using the field OBJECTID (standard in fgdb fc's) and the field is not delimited. To make this a little more gentle the code should be changed a bit: 1) use the describe methode to obtain the oid field (if it has one, if not exit the procedure): desc = arcpy.Describe(fcpoints) if desc.hasOID:     fldOID = desc.OIDFieldName else:     exit("featureclass does not have an Object ID field") 2) the whereclause should be build differently:         fld = arcpy.AddFieldDelimiters(fcpoints,fldOID)         whereclause = "{1} = {0}".format(oid,fld) This will result probably in a whereclause like: "FID" = 0 So the final code in your case would be something like: import arcpy, os, numpy from numpy import * # settings of your data fcpoints = r'G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\clip16points.shp' fldVisible = 'VisiblePix' # will contain the number of visible pixels dem = r'G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\combraster2' # ouput workspace and name for viewshed outputWS = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb' prefixVS = 'Viewshed' arcpy.env.workspace = outputWS arcpy.env.overwriteOutput = True # grab a SA license if arcpy.CheckExtension("Spatial") == "Available":     arcpy.CheckOutExtension("Spatial") else:     exit("no SA license available") # get the oid field desc = arcpy.Describe(fcpoints) if desc.hasOID:     fldOID = desc.OIDFieldName else:     exit("featureclass does not have an Object ID field") # loop through points fields = ("oid@",fldVisible) with arcpy.da.UpdateCursor(fcpoints, fields) as cursor:     for row in cursor:         oid = row[0]         # define output name (of oid = 1, raster will be called Viewshed1         outVSname = outputWS + os.sep + prefixVS + str(oid)         # create a whereclause         fld = arcpy.AddFieldDelimiters(fcpoints,fldOID)         whereclause = "{1} = {0}".format(oid,fld)         # create featurelayer for current observer point         FLsources = "sources_Layer"         arcpy.MakeFeatureLayer_management(fcpoints,FLsources,whereclause,outputWS)         # calculate viewshed         arcpy.gp.Viewshed_sa(dem,FLsources,outVSname,"1","FLAT_EARTH","0.13")         # determine number of visible cells         myArray = arcpy.RasterToNumPyArray(outVSname)         myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape))         VisiblePix = myArray2.sum()         print "There are {0} visible pixels for point with oid={1}".format(VisiblePix,oid)         # update value in point featureclass         row[1] = VisiblePix         cursor.updateRow(row)         del myArray, myArray2 del row # release license arcpy.CheckInExtension("Spatial") Before you run it please check you input dem and environment settings as explained in my previous post. It this still doesn't work, try to run the code with a local copy of your dem. The 'G' drive most likely maps to a network folder and I have seen issues with running scripts with grids on network locations that failed, when running it with a local copy actually worked. Kind regards, Xander ... View more

Hi Anne, The description of this error is: All of the cells in the grid are NoData. The tool requires that some cells have values other than NoData to operate on. This may happen if there is no overlap between the extent and the input(s). Solution Ensure that the input grid has some cell values that are not NoData. If you have set the Extent environment, make sure it falls at least partway within your input data. Can you verify that this is not the case? Kind regards, Xander ... View more

To lend a hand: Perform a spatial join from your points to your polygons I assume your polygon fc has a numeric field which is unique for each polygon Define this field as the fldID This code will loop over each polygon, read out the geometry and set up a where clause expression. This is fed to the disperse function where only the points in the polygon are dispersed. This is a result of the code: [ATTACH=CONFIG]29470[/ATTACH] Another option is to change the disperse function and initially determine if the punt is within the current polygon and process the point only if it is inside (don't affect the point location if it's not within the current polygon)... #------------------------------------------------------------------------------- # Name:        Disperse2.py # Purpose:     Disperse points in multiple polygons # Author:      arcpy Team #              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/ # Created:     29-11-2013 #------------------------------------------------------------------------------- def main():     global arcpy, random     import arcpy, random     fcPoints = r"C:\Project\_Forums\Disperse\test.gdb\points_joined" # points after spatial join     fcPolygon = r"C:\Project\_Forums\Disperse\test.gdb\Polygons"     fldID = 'myID' # polygon ID field in point FC after spatial join     arcpy.env.overwriteOutput = True     with arcpy.da.SearchCursor(fcPolygon, ("SHAPE@",fldID)) as cursor:         for row in cursor:             polygon = row[0]             myId = row[1]             expression = arcpy.AddFieldDelimiters(fcPoints, fldID) + " = {0}".format(myId)             disperse_points(fcPoints,polygon,expression)     del row     print "ready..." def point_in_poly(poly, x, y):     pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)     return poly.contains(pg) def disperse_points(in_points, polygon, expression):     lenx = polygon.extent.width     leny = polygon.extent.height     with arcpy.da.UpdateCursor(in_points, "SHAPE@XY", where_clause=expression) as points:         for p in points:             x = (random.random() * lenx) + polygon.extent.XMin             y = (random.random() * leny) + polygon.extent.YMin             inside = point_in_poly(polygon, x, y)             while not inside:                 x = (random.random() * lenx) + polygon.extent.XMin                 y = (random.random() * leny) + polygon.extent.YMin                 inside = point_in_poly(polygon, x, y)             points.updateRow([(x, y)]) if __name__ == '__main__':     main() Kind regards, Xander ... View more

Hi Peter, I changed my code from this thread a bit and this should create the rectangles: import arcpy # input point file fc = r'C:\Project\_Forums\pointInterpol\fgdb\test.gdb\points' # edit this fields = ("SHAPE@XY") # size of rectangle h = 1000 # height, edit this w = 1500 # width, edit this # output feature class fcout = r'C:\Project\_Forums\pointInterpol\fgdb\test.gdb\diamonds01' # edit this features = [] with arcpy.da.SearchCursor(fc, fields) as cursor:     for row in cursor:         arrPnts = arcpy.Array()         X = row[0][0]         Y = row[0][1]         # point 1         pnt = arcpy.Point(X-w/2,Y-h/2)         arrPnts.add(pnt)         # point 2         pnt = arcpy.Point(X-w/2,Y+h/2)         arrPnts.add(pnt)         # point 3         pnt = arcpy.Point(X+w/2,Y+h/2)         arrPnts.add(pnt)         # point 4         pnt = arcpy.Point(X+w/2,Y-h/2)         arrPnts.add(pnt)         # point 5 (close diamond)         pnt = arcpy.Point(X-w/2,Y-h/2)         arrPnts.add(pnt)         pol = arcpy.Polygon(arrPnts)         features.append(pol) # write to output arcpy.CopyFeatures_management(features, fcout) Kind regards, Xander ... View more

Hello. I was wondering if I could get a bit of help with this code? I've been trying to run it in PyScripter. I'm getting the following error message:         TypeError: unsupported operand type(s) for ^: 'int' and 'numpy.float64' Not sure if this helps, but in this thread: http://forums.arcgis.com/threads/54169-PyScripter-Question ... they say that you need to use the 32 bits version of PyScripter. Kind regards, Xander ... View more

The points getting placed far away from any ZIP code polygon may be due to differences in coordinate system. Make sure your points and ZIP code polygons share the exact same coordinate system. The point coordinates are calculated using the extent of a single polygon. Using multiple polygons would require some changes. You would have to loop through all the polygons and create a feature layer based on a select by location in which only those points are passed through that are within the current polygon. Another option would be a spatial join from the points to the polygons. By this you will probably obtain the FIPS at the points. Use this FIPS in a whereclause to obtain the right polygon and also in the whereclause to process only those points within that polygon. Kind regards, Xander ... View more

I tested with this code and it works as I would expect. The point is that you need to provide a polygon geometry and not a polygon featureclass to the disperse_points function. #------------------------------------------------------------------------------- # Name:        Disperse.py # Purpose:     Disperse points in polygon # Author:      arcpy Team #              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/ # Created:     27-11-2013 #------------------------------------------------------------------------------- def main():     global arcpy, random     import arcpy, random     fcPoints = r"C:\Project\_Forums\Disperse\test.gdb\points" # edit this     fcPolygon = r"C:\Project\_Forums\Disperse\test.gdb\polygon" # edit this     with arcpy.da.SearchCursor(fcPolygon, "SHAPE@") as cursor:         for row in cursor:             polygon = row[0] # take first polygon             break     del row     disperse_points(fcPoints,polygon)     print "ready..." def point_in_poly(poly, x, y):     """Returns if the point is inside the polygon.     Parameters:         poly: arcpy.Polygon() geometry         x:    x coordinate (float)         y:    y coordinate (float)     """     pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)     return poly.contains(pg) def disperse_points(in_points, polygon):     """Randomly disperse points inside a polygon.     Parameters:         in_points:  Point feature class/layer (with or without selection)         polygon:    arcpy.Polygon() geometry     """     lenx = polygon.extent.width     leny = polygon.extent.height     with arcpy.da.UpdateCursor(in_points, "SHAPE@XY") as points:         for p in points:             x = (random.random() * lenx) + polygon.extent.XMin             y = (random.random() * leny) + polygon.extent.YMin             inside = point_in_poly(polygon, x, y)             while not inside:                 x = (random.random() * lenx) + polygon.extent.XMin                 y = (random.random() * leny) + polygon.extent.YMin                 inside = point_in_poly(polygon, x, y)             points.updateRow([(x, y)]) if __name__ == '__main__':     main() Kind regards, Xander ... View more

Dear Xander, Thank you so much for the code ! It worked in perfect manner :). Hi Milan, I'm glad I could be of any assistance. If you think it was helpful you    can use the "arrow" button in order to help other  members find useful   information: More info here: http://resources.arcgis.com/en/help/forums-mvp/ Kind regards, Xander ... View more

Hi Anne, I think that you can use python for this. You would probably have to convert the model to python and clean it up a bit. I suppose you loop through the observer points and for each point you are calculating a viewshed raster. This raster contains 1 for visible and 0 for not visible. Using numpy it is easy to get the sum of the raster. In case you have NoData cells, the numpy array needs to be correct for this. The code is very simple: import numpy, arcpy from numpy import * myArray = arcpy.RasterToNumPyArray("Viewshed1") # avoid erroneous sum due to NoData cells myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape)) visibleCells = myArray2.sum() In the output field of the point featureclass (observers), you can use a simple updatecursor (and a whereclause) to set the calculated value. The script could look something like this: import arcpy, os, numpy from numpy import * # settings of your data fcpoints = r'C:\Project\_Forums\Viewshed\test.gdb\sources' fldVisible = 'VisiblePix' # will contain the number of visible pixels dem = r'C:\Project\_Forums\cost\fgdb\test.gdb\elevutm16n' # ouput workspace and name for viewshed outputWS = r'C:\Project\_Forums\Viewshed\test.gdb' prefixVS = 'Viewshed' arcpy.env.workspace = outputWS arcpy.env.overwriteOutput = True # grab a SA license if arcpy.CheckExtension("Spatial") == "Available":     arcpy.CheckOutExtension("Spatial") else:     exit("no SA license available") # loop through points fields = ("OID@",fldVisible) with arcpy.da.UpdateCursor(fcpoints, fields) as cursor:     for row in cursor:         oid = row[0]         # define output name (of oid = 1, raster will be called Viewshed1         outVSname = outputWS + os.sep + prefixVS + str(oid)         # create featurelayer for current observer point         whereclause = "OBJECTID = {0}".format(oid)         FLsources = "sources_Layer"         arcpy.MakeFeatureLayer_management(fcpoints,FLsources,whereclause,outputWS)         # calculate viewshed         arcpy.gp.Viewshed_sa(dem,FLsources,outVSname,"1","FLAT_EARTH","0,13")         # determine number of visible cells         myArray = arcpy.RasterToNumPyArray(outVSname)         myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape))         visPixels = myArray2.sum()         print "There are {0} visible pixels for point with OID={1}".format(visPixels,oid)         # update value in point featureclass         row[1] = visPixels         cursor.updateRow(row)         del myArray, myArray2 del row # release license arcpy.CheckInExtension("Spatial") If you have any questions, please let me know. Kind regards, Xander ... View more