Latest Contributions by ChrisSnyder

How about the euclidean alocation tool? (Spatial Analyst extension required) http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z0000001m000000.htm Maybe some form of the buffer tool output intersected with vectorized euclidean alocation output? ... View more

You said you are running this on a Linux box... I don't have any experience with Server runninng on Linux, but can you start up a Python window (run Python.exe) and run the command: import arcgisscripting or import arcpy ? If not, this either indicates: 1. There is something wrong with your ArcServer/Python install 2. You can't run geoprocessing (gp/arcpy) scripts via ArcServer on Linux (I think you can though if I'm not mistaken??) The error messages you get seem to hint at the 1st possibility. ... View more

So if I am not mistaken, you have a DEM with some small holes in it... Maybe LiDAR data from some fly-by-night vendor? And you want to fill each hole with "reasonable" elevation values. I think what you actually might want to do is for each little nodata hole, get the minimum elevation value (maybe the MEAN?) of a 1 cell buffer zone around each hole... each hole probably has a different elevation values surounding it, so you want to be sure to use the RegionGroup tool as well. This expresion seemed to work for me... Just replace "my_grid" with the name of your DEM. It will produce a "no hole" DEM where each hole has the minimum elevation value of 1 pixel buffer surounding the hole. noHolesGrd = Con(IsNull("my_grid"), ZonalStatistics(FocalStatistics(RegionGroup(Con(IsNull("my_grid"), 1)), NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA"), "VALUE", "my_grid", "MINIMUM", "DATA"), "my_grid") Be aware that this expresion will also "fill in" no data areas on the periphery of your dataset as well, and that might not be desirable... The code to exclude these periphery areas would be more complicated... Otherwise you could also manually create a mask layer specifically showing the areas you intend to fill (edit the mask polygon so as to delete the periphery nodata areas). ... View more

Are you sure the FGDB isn't compressed (aka in a read only state)? How large is a "larger dataset"? ... View more

If anyone uses any of that code above, I found a bit of a bug.. Basically if your child subprocess is quite chatty (like a bunch of logging to a text file), the subprocess.stdout buffer in the memory can fill up and cause the child subprocess to hang forever! If stdout is sent to PIPE and is more than 65536 characters, the subprocess will hang. After pulling my hair out for a few hours, I found the answer here: http://thraxil.org/users/anders/posts/2008/03/13/Subprocess-Hanging-PIPE-is-your-enemy/ Some work arounds: -------------------- 1. Don't have a chatty subproces script (stdout memory buffer is 65536 characters or so, which can fill up fast). 2. Don't use stdout = PIPE if you don't care about stdout messages. This is what I did in my case, because all my parent script needs is the subprocess .poll() and .returncode So for example, I did this: subprocess.Popen([jobDict[jobId][0], jobDict[jobId][1], str(inputVar1), str(inputVar2)], shell=False) instead of this: subprocess.Popen([jobDict[jobId][0], jobDict[jobId][1], str(inputVar1), str(inputVar2)], shell=False, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE) 3. Send stdout to an actual file on disk instead of memory... all the "filelike" stdout methods (like .readlines) are the same anyway! File or memory - doesn't really matter. ... View more

No need to loop - just this: import arcpy, time arcpy.env.workspace = r'C:\TempArcGISCalculate\TempFile\TestforReadingFile' arcpy.env.overwriteOutput=True ptfc = 'Section01_SplitedTrips.shp' t0=time.clock() bufferFC = 'in_memory/Section01Buffer' arcpy.Buffer_analysis(ptfc, bufferFC,'4 meters') identityFC = r"in_memory\identity" arcpy.Identity_analysis(ptfc, bufferFC, identityFC, "ALL", "", "") Then I think all the info you should need is it the attribute table of  identityFC Maybe add/populate an X andY field. Beware that you will now have poiuts that will overlap n times in areas that have n overlapping buffer polygons. Then use a searchcursor and a python dictionary. You can key the dictionary by (x,y) coordinate pairs. Might have to build a few different dictionaries to keep track of everything. ... View more

Thanks for the narative. So I think you can: 1. Buffer your points 2. Run an Identity (or Intersect) tool on your points and buffer polygon - This will "tag" each point with what point buffer it is within - be carefull because you will end up with duplicate points where the buffer polygons overlap!!! 3. In a cursor (maybe helpfull to use a Python dictionary to sort out some of the realtionships) - use the Intersect/Identity tabular information (x/y/point_id/buffer_id) to analyze the realtionships. Run a small sample area of your point and point buffers through the Identity or Intersect tools to see what can be done with the output table info. A recursive SelectByLocation will be extreemly slow compared with the convetional Intersect or Identity tools. ... View more

And you need to import the time module to get time.clock() to work import time ... View more

Anything is possible... Can you give a narative description of what this code is doing? Looks like it is attempting to define some sort of 'spatial association' - like identifying the points within the buffer of a polygon. Maybe a conventional overlay tool would better serve this purpose... For example, Intersect (Intersect tool) your polygons buffers with your points and cursor through that output. ... View more

or better yet... exampleDict = dict([(r.KEY, (r.VALUE1, r.VALUE2)) for r in arcpy.SearchCursor(myTbl)]) ... View more

Right on... I think I'm getting into the groove now. Thanks a million Phil! ... View more

New strugle: I had written some code years ago that would create a 1 cell outline of the DATA portion of a raster (see attached image): somaExp = "con(" + flowDirGrd + " > 0 & isnull(Focalmax(" + flowDirGrd + ", 'RECTANGLE', 3, 3, 'NODATA')), 1, setnull(" + flowDirGrd + "))" ....where 'flowDirGrd' is just some raster (doesn't neccesarily have to be a flow direction raster). So basically just the outer DATA cells get coded as a 1 and everything else is NoData. I came up with this for v10 sytax: demRingTmp = arcpy.sa.Con(flowDirTmp > 0 & arcpy.sa.IsNull(arcpy.sa.FocalStatistics(flowDirTmp, arcpy.sa.NbrRectangle(3, 3, "CELL"), "MAXIMUM", "NODATA")), 1) ...which runs but produces the wrong results for some reason. Anyone wanna take a stab at writting this ugly thing in v10 syntax (and get the correct result)? [ATTACH=CONFIG]12500[/ATTACH] ... View more

This code should give you a random 10% sample of the OBJECTIDs in a feature layer that have a selected set. #Warning: Untested import random, arcpy samplePct = 0.1 arcpy.MakeFeatureLayer_managment(roadsFC, "roads", "") arcpy.SelectLayerByAttribute_managment("roads", "NEW_SELECTION", "ROAD_TYPE = 'Highway') fidSetList = [int(fid) for fid in arcpy.Describe("roads").fidSet.split(";")] if len(fidSetList) * samplePct >= 1: #error check rto make sure you have enough selected features to get an appropriate sample    randomSampleList = radom.sample(fidSetList, int(len(fidSetList) *samplePct)) else:    print "Not enough selected features to constitute a " + str(samplePct * 100) +  "% sample!" Then you could delete them like this: arcpy.MakeFeatureLayer(roadsFC, "delete_me", "OBJECTID in (" + str(randomSampleList)[1:-1] + ")") arcpy.DeleteFeatures_Managment("delete_me") or use an update cursor... same thing. ... View more

Thanks Phil... I was missing the whole "casting as a raster" thing... and the implied setnull - didn't know about that - very usefull! Took me a while to realize that the temporary raster objects (for example "slopePctTmp" in the expresion: slopePctTmp = arcpy.sa.Slope(fillGrd, "PERCENT_RISE", "") can and should be used even after you ".save()" them, and low and behold, they now reference the updated .save() path! One you you use a variable to reference the .save() output path, that variable is no longer a "raster object", so just keep using the raster reference (slopePctTemp). The big lightbulb though - and a **MAJOR ARCPY CRASH BUG ISSUE***: What I had been doing: Seems that you CAN use the the output path variable (as opposed to the raster object variable) in SINGLE tool expreessions... for example: #Process: Run a Fill fillTmp = arcpy.sa.Fill(clipDem, "") fillGrd = areaIdFolderPath + "\\fill" fillTmp.save(fillGrd) #Process: FlowDir flowDirTmp = arcpy.sa.FlowDirection(fillGrd, "", "") #NOTE: I am using the fillGrd variable not the the fillTmp variable: BOTH work for "single tool expressions" flowDirGrd = areaIdFolderPath + "\\flowdir" flowDirTmp.save(flowDirGrd) But if you try to use the output path variables (and not the raster object variables) in MULTI TOOL EXPRESSIONS Python will CRASH BOOM BLOWUP... for example: This is OK: #Process: FlowAcc flowAccTmp = arcpy.sa.FlowAccumulation(flowDirGrd, "", "INTEGER"); showGpMessage() flowAccGrd = areaIdFolderPath + "\\flowacc" flowAccTmp.save(flowAccGrd) #Process: Create a streamlink acreThreshold = 1 acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2) flowAccGrd = arcpy.Raster(flowAccGrd) streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccTmp > acreThresholdInPixels, 1), flowDirGrd) streamLinkGrd = areaIdFolderPath + "\\strmlnk" streamLinkTmp.save(streamLinkGrd) This VERY bad and cause Python "Application Error" explosion: #Process: FlowAcc flowAccTmp = arcpy.sa.FlowAccumulation(flowDirGrd, "", "INTEGER"); showGpMessage() flowAccGrd = areaIdFolderPath + "\\flowacc" flowAccTmp.save(flowAccGrd) #Process: Create a streamlink acreThreshold = 1 acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2) flowAccGrd = arcpy.Raster(flowAccGrd) streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccGrd > acreThresholdInPixels, 1), flowDirGrd) streamLinkGrd = areaIdFolderPath + "\\strmlnk" streamLinkTmp.save(streamLinkGrd) So the best solution is to ***ALWAYS*** use the raster objects andnot the path variables: #Process: FlowAcc flowAccTmp = arcpy.sa.FlowAccumulation(flowDirGrd, "", "INTEGER"); showGpMessage() flowAccGrd = areaIdFolderPath + "\\flowacc" flowAccTmp.save(flowAccGrd) #Process: Create a streamlink acreThreshold = 1 acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2) flowAccGrd = arcpy.Raster(flowAccGrd) streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccTmp > acreThresholdInPixels, 1), flowDirTmp) streamLinkGrd = areaIdFolderPath + "\\strmlnk" streamLinkTmp.save(streamLinkGrd) So, moral of the story... Don't use the .save(pathvariable) path variables in raster expresions!!! ... View more