|
POST
|
I got this to work for at least tifs and shapefiles. I'm not sure why I had to use mapping.layer, line 20 and 29, before I could get an Extent of the data. Seems like the overlap tests, lines 32 - 47, could be down more efficiently. Is there just one test that would find if there is no touching or overlap, then skip that file. Where's the docs on that .overlaps stuff? Can't find it today. Currently all my output lands in one folder. My next step needs to be writing the output to separate folders in the output directory, making the output tree structure look like the input tree. I'm guessing that will involve manipulating the strings of the full paths of the inputs.
#testclip.py
#test trying to walk a directory tree, find rasters and shapefiles, and clip to a area of Interest (AOI)
# and make it run in 10.0 so without da.walk.
# Paul Huffman, Yakama Fisheries, 7/7/2014
# Import arcpy module
import arcpy
import os
#arcpy.CheckOutExtension("spatial")
# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "c:\\avdata\\PythonTest\\ClipLidar\\Scratch"
arcpy.env.workspace = "c:\\avdata\\PythonTest\\ClipLidar\TestData"
arcpy.env.overwriteOutput = True
clipshape = "c:\\avdata\\PythonTest\\ClipLidar\\AOIPoly_SP.shp"
InputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input"
OutputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Output"
#Find AOI extent
AOI = arcpy.mapping.Layer(clipshape)
AOIextent=AOI.getExtent()
#Find the input feature classes
for root, dirs, files in os.walk(InputFolder):
for name in files:
Subset = 0
if os.path.splitext(name)[1] == ".tif":
#print os.path.join(root, name), '\n'
FileObj = arcpy.mapping.Layer(os.path.join(root, name))
FileExtent = FileObj.getExtent()
print "Processing: " + name
Extent_Overlap = str(FileExtent.overlaps(AOIextent))
Extent_Within = str(FileExtent.within(AOIextent))
Extent_Touches = str(FileExtent.touches(AOIextent))
Extent_Crosses = str(FileExtent.crosses(AOIextent))
if Extent_Overlap == 'True':
Subset = 1
print "Extent_Overlaps"
elif Extent_Within == 'True':
Subset = 1
print "Extent_Within"
elif Extent_Touches == 'True':
Subset = 1
print "Extent_Touches"
elif Extent_Crosses == 'True':
Subset = 1
print "Extent_Crosses"
if Subset == 1:
arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputFolder,name)),clipshape,"#","ClippingGeometry")
if os.path.splitext(name)[1] == ".shp":
print(os.path.join(root, name))
arcpy.Clip_analysis((os.path.join(root, name)),clipshape,(os.path.join(OutputFolder,name)))
... View more
07-08-2014
02:59 PM
|
0
|
10
|
2604
|
|
POST
|
I was trying to get this script done yesterday, but puzzled why half of it doesn't work. Because ESRI is doing site maintenance, I couldn'tget back in to ask this question. Client still waiting. I was trying to get a script to walk down a tree of LIDAR data, find any rasters and feature classes, and clip them to an area of interest. My test data tree just has some tifs in one folder, some shapefiles in another to make things simpler. And I need to make this run in ArcGIS 10.0, so da.walk and da.Cursor are not available. When my script hit the first tif file, the AOI polygon didn't overlap it at all and I got an error message that it didn't overlap. (Failed to create raster dataset. The reason could be: The clip feature is outside the raster extent. Failed to execute (Clip).) I was going to put off solving the problem of testing for overlap until I could get the file finding and clipping figured out. So, following an example at http://gis.stackexchange.com/questions/13571/arcpy-compare-extents-of-rasters, I added the commented out section, lines 26 to 49, to check for overlap of extent for the tifs. However, the script hits line 48, and says name 'Subset' is not defined. Also none of the print statements in the section are executed, so it looks like none of the extent tests worked. What went wrong? If I just run the script down to the first print statement, all the tifs are found and printed. In the next section, an if statement looks for just shapefiles, and successfully clips them and writes the output to a new output tree. No check for overlap is done for the shapefiles but I probably should to eliminate empty output shapefile and to save execution time. Hey, how do you insert code in this new forum? #testclip.py #test trying to walk a directory tree, find rasters and shapefiles, and clip to a area of Interest (AOI) # and make it run in 10.0 so without da.walk. # Paul Huffman, Yakama Fisheries, 7/7/2014 # Import arcpy module import arcpy import os arcpy.CheckOutExtension("spatial") # Set Geoprocessing environments arcpy.env.scratchWorkspace = "c:\\avdata\\PythonTest\\ClipLidar\\Scratch" arcpy.env.workspace = "c:\\avdata\\PythonTest\\ClipLidar\TestData" arcpy.env.overwiteOutput = True clipshape = "c:\\avdata\\PythonTest\\ClipLidar\\AOIPoly.shp" InputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input" OutputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Output" #Find AOI extent AOI = arcpy.mapping.Layer(clipshape) AOIextent=AOI.getExtent() #Find the input feature classes for root, dirs, files in os.walk(InputFolder): for name in files: if os.path.splitext(name)[1] == ".tif": print os.path.join(root, name), '\n' FileObj = arcpy.mapping.Layer(os.path.join(root, name)) FileExtent = FileObj.getExtent() print "Processing: " + name print FileExtent ## Extent_Overlap = str(FileExtent.overlaps(AOIextent)) ## Extent_Within = str(FileExtent.within(AOIextent)) ## Extent_Touches = str(FileExtent.touches(AOIextent)) ## Extent_Crosses = str(FileExtent.crosses(AOIextent)) ## if Extent_Overlap == 'True': ## Subset = 1 ## print "Extent_Overlaps" ## elif Extent_Within == 'True': ## Subset = 1 ## print "Extent_Within" ## elif Extent_Touches == 'True': ## Subset = 1 ## print "Extent_Touches" ## elif Extent_Crosses == 'True': ## Subset = 1 ## print "Extent_Crosses" ## #print(os.path.join(root, name)) ## #arcpy.Clip_management(name,clipshape,(os.path.join(OutputFolder,name))) ## if Subset == 1: ## arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputFolder,name)),clipshape,"#","ClippingGeometry") ## if os.path.splitext(name)[1] == ".shp": ## print(os.path.join(root, name)) ## #arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputFolder,name)),clipshape,"#","ClippingGeometry") ## arcpy.Clip_analysis((os.path.join(root, name)),clipshape,(os.path.join(OutputFolder,name))) End of code. I also found the python glob method. It also can find all the tifs by searching through the tree like: #Find the input feature classes print (glob.glob("c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input\\*\\*.tif"),'\n') ##for root, dirs, files in os.walk(InputFolder): ## for name in files: ## if os.path.splitext(name)[1] == ".tif": ## print os.path.join(root, name), '\n' where the glob statement returns the same list as the next four commented out statements. However, I eventually want to parse out the folders in which the files were found so I can reconstruct the input tree structure in the output tree. Don't know I can do that with the glob returns. Actually I don't know how to to that yet with the walk returns either.
... View more
07-08-2014
10:06 AM
|
0
|
11
|
2604
|
|
POST
|
It looks like I need to do it in 10.0 unless I travel there with a > 10.0 laptop installation. And I might have to do that because some of that data is in geodatabases, so os.walk will not find them. Something else has come up. Many of the input feature classes are way outside the AOI clip. Too much time is being taken up by handling these out of range FCs and the output files are just empty. Is there a way of avoiding the out of range data in the clipping loop? It looks like I could first test for intersection with the clip polygon with Select layer by Location but I'm having a hard time telling if that tool is available at 10.0. I don't have 10.0 installed anywhere at my current location, and I was unable to find 10.0 help online.
... View more
07-02-2014
12:42 PM
|
0
|
0
|
2604
|
|
POST
|
I want to loop through a folder and clip all feature classes by a AOI poly and write the results to an output folder. How to I also have the loop search all subfolders for vectors and raster feature classes as well?
... View more
07-02-2014
09:01 AM
|
0
|
19
|
10580
|
|
POST
|
Lately I've been trying to figure out how much longer my storage on my ArcServer is going to last. The major portion of the data storage is LIDAR data and its products that we procure every year. I thought maybe one approach would be to see how much each year's LIDAR set is using, then I might be able to guess how many more years of LIDAR data I can fit on the remaining disk space. But because the data is all packed away in SDE storage, I'm having difficulty figuring out how I can read disk usage for each year. The only way I can figure out how to do it is backup the database, compact, look at remaining free disk space, delete a year's set, compact, then look at the increase in free disk space, then restore. Are there any other good ideas? Arc Server 9.3.1 on a Windows Server platform, MS SQL Server
... View more
06-25-2014
02:50 PM
|
0
|
1
|
752
|
|
POST
|
nidhinkn - why is that command line step to copy connection file from toolbox to Catalog needed?
... View more
06-03-2014
02:20 PM
|
0
|
0
|
856
|
|
POST
|
ESRI tech support did a screen share with me and quickly found the problem. The "enable background processing" check box was on in the 2012 map document. Even though I could "print env.workspace" and get the correct path, background processing somehow resets everything, apparently, on execution of a next line so the workspace path and the input shapefile can't be found. ESRI thought they might submit a bug report on it, since this doesn't seem to be desired or expected behavior. From all the posts in this forum, it seems like "enable background processing" is often involved in unexpected behavior in many different situations.
... View more
05-30-2014
01:43 PM
|
0
|
0
|
1109
|
|
POST
|
The guy I was helping with this told me that he had tried using the field calculator, but was overwhelmed by its complexity, so I steered away from the field calculator too. Now he has two ways to do this. Using da.cursor without opening ArcMap will be nice with large data sets, but I couldn't notice the difference with my little test data set. The longest part seemed to taken just loading arcpy. My workflows often involve having ArcMap open so I can view and edit feature classes as I run through them. But that doesn't mean I can't drag some da.cursor code into the ArcMap python window for this kind of calculation.
... View more
05-28-2014
01:23 PM
|
0
|
0
|
2083
|
|
POST
|
What's wrong with this approach, using Field Calculator with a simple python function, no cursors? Works well on my small test set. [ATTACH=CONFIG]34118[/ATTACH]
... View more
05-27-2014
02:07 PM
|
0
|
0
|
2083
|
|
POST
|
Thanks guys. All I needed to add was del row
del rows to release the locks on the attribute table.
... View more
05-27-2014
09:58 AM
|
0
|
0
|
2083
|
|
POST
|
Apparently I just need to add a line cursor.updateRow(row) at the end of the loop to write the changes. (As shown in the documentation, if I would ever just read farther along. ) But another problem showed up. I was trying to show a colleague how to do this, but he only has ArcGIS 10.0, so he doesn't have da.UpdateCursor. I need to come up with a way using just arcpy.UpdateCursor. I've tried a couple different loops but come up with either the input file doesn't exist or the field field1 doesn't exist.
# Import arcpy module
import arcpy
from arcpy import env
env.workspace = "C:\\avdata\\PythonTest\\testgdb.gbd\\"
# Local variables:
fc = "zilltogra" #input Feature class
field1 = "PassNum" #field with test values
field2 = "ReddStatus" #field to be calculated
items = arcpy.ListFields(fc)
print items
##cursor = arcpy.UpdateCursor(inshape)
##row = cursor.next()
##while row:
## if row.field1 == 1:
## row.field2 = "Definite"
## elif row.field1 == 2:
## row.field2 = "Probable"
## cursor.updaterow(row)
## row = cursor.next(1)
rows = arcpy.UpdateCursor(fc)
for row in rows:
if row.field1 == 1:
row.field2 = "Definite"
elif row.field1 == 2:
row.field2 = "Probable" # Import arcpy module
import arcpy
from arcpy import env
env.workspace = "C:\\avdata\\PythonTest\\"
# Local variables:
inshape = "zilltogra101013.shp" #input shapefile
field1 = "PassNum" #field with test values
field2 = "ReddStatus" #field to be calculated
fields = ['PassNum', 'ReddStatus']
##cursor = arcpy.da.UpdateCursor(inshape)
##row = cursor.next()
##while row:
## if row.field1 == 1:
## row.field2 = "Definite"
## elif row.field1 == 2:
## row.field2 = "Probable"
## cursor.updaterow(row)
## row = cursor.next(1)
rows = arcpy.UpdateCursor(inshape)
for row in rows:
if row.field1 == 1:
row.field2 = "Definite"
elif row.field1 == 2:
row.field2 = "Probable"
... View more
05-23-2014
03:35 PM
|
0
|
0
|
2083
|
|
POST
|
I was testing out arcpy.da.UpdateCursor to calculate a field based on another fields value. My script runs without error, but doesn't seem to make any changes in the shapefile. I tried adding the line overwriteOutput = "True" but that didn't work. The print line I added for debugging shows the correct values for ReddStatus or row[1], but they don't get saved in the output. # Import arcpy module import arcpy from arcpy import env env.workspace = "C:\\avdata\\PythonTest\\" arcpy.env.overwriteOutput = "True" # Local variables: inshape = "zilltogra101013.shp" #input shapefile field1 = "PassNum" #field with test values field2 = "ReddStatus" #field to be calculated fields = ['PassNum', 'ReddStatus'] with arcpy.da.UpdateCursor(inshape, fields) as cursor: for row in cursor: if row[0] == 1: row[1] = "Definite" elif row[0] == 2: row[1] = "Probable" print row[0],row[1], '\n'
... View more
05-23-2014
03:22 PM
|
0
|
7
|
2296
|
|
POST
|
What makes this more weird is I was running a script much like this two weeks ago in the python window in a different map document. I opened up the 2013 data map document again and tried it. Yes, this script still runs, finds the shapefile just fine, just a different workspace, different map, same terrible long shapefile names. There must be something different about the 2012 map. Can there be something wrong with the workspace?
# Import arcpy module
import arcpy
from arcpy import env
env.workspace = "C:\\avdata\\FallChinookRedds\\2013\\Shapefiles\\"
# Local variables:
inshape = "MDhighway972211113.shp"
fname = "\""+inshape+"\""
surdate = "\"11/01/2013\""
streamname = "\"Yakima River\""
specname = "\"Fall Chinook\""
# Process: Delete Field
arcpy.DeleteField_management(inshape, "TYPE;Y_PROJ;X_PROJ;DISPLAY;SYMBOL;UNUSED1;DIST;PROX_INDEX;COLOR;DEPTH;TEMP;TIME;WPT_CLASS;SUB_CLASS;ATTRIB;LINK;STATE;COUNTRY;CITY;ADDRESS;FACILITY;CROSSROAD;UNUSED2;ETE;DTYPE;MODEL")
# Process: Add Field
#arcpy.AddField_management(inshape, "COMMENT", "TEXT", "", "", "50", "", "NON_NULLABLE", "NON_REQUIRED", "")
#arcpy.AddField_management(inshape, "ALTITUDE", "DOUBLE", "", "", "8", "", "NON_NULLABLE", "NON_REQUIRED", "")
#arcpy.AddField_management(inshape, "FILENAME", "TEXT", "", "", "254", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "SurveyDate", "DATE", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "SurveyYear", "LONG", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "observer", "TEXT", "", "", "25", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "PassNum", "SHORT", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "Species", "TEXT", "", "", "50", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "Stream", "TEXT", "", "", "50", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(inshape, "ReddStatus", "TEXT", "", "", "3", "", "NON_NULLABLE", "NON_REQUIRED", "")
# Process: Calculate Field
arcpy.CalculateField_management(inshape, "SurveyDate", surdate, "VB", "")
arcpy.CalculateField_management(inshape, "SurveyYear", "2013", "VB", "")
arcpy.CalculateField_management(inshape, "Stream", streamname, "VB", "")
#arcpy.CalculateField_management(inshape, "observer", "\"Brandon Rogers\"", "VB", "")
#arcpy.CalculateField_management(inshape, "PassNum", "PASS", "VB", "")
arcpy.CalculateField_management(inshape, "Species", specname, "VB", "")
arcpy.CalculateField_management(inshape, "FILENAME", fname, "VB", "")
arcpy.GetCount_management(inshape)
... View more
05-07-2014
01:49 PM
|
0
|
0
|
1109
|
|
POST
|
Since you are playing the dangerous game of working parallel with two apps in the same workspace, I'm thinking you may want to try this method to have arcpy verify what's REALLY there in the workspace, just not last time it checked: arcpy.RefreshCatalog(env.workspace) Didn't help.
... View more
05-06-2014
09:51 AM
|
0
|
0
|
1109
|
|
POST
|
Today I've been successfully cranking through these shapefiles with this script. Making so much progress I started to think that yesterday I just had a typo in the workspace path that I didn't catch. But now I see that it doesn't matter if the shapefile name is 13 characters long. It doesn't matter if the shapefile name starts with a numeral. What causes my script to fail is running it in the ArcMap Python window. It works fine in the ArcCatalog python window, or from IDLE if I shut down ArcMap and ArcCatalog. In ArcMap, is there some env variable that persists even after the script sets the env.workspace variable? Does it revert to the map properties right after the script runs the env.workspace line? update: I looked again at the map document properties. Home folder was set correctly, the default geodatabase was set to a different location, so I made a test file geodatabase in the home directory, set it as default. Then I examined the map document geoprocessing options, and made sure that the workspace was set to the home folder. The home folder contains all the shapefiles and a geodatabase is not used in the map. Yet when I try to execute any geoproccessing line in my code in the ArcMap python window, I get "ERROR 000732: Input Table: Dataset MarionHarrahTOoldgold102512.shp does not exist or is not supported Failed to execute (CalculateField)." as if the path is wrong. This is code that works in ArcCatalog or IDLE. I don't know what is different about ArcMap python window.
... View more
05-05-2014
03:16 PM
|
0
|
0
|
1109
|
| Title | Kudos | Posted |
|---|---|---|
| 1 | 11-28-2018 04:57 PM | |
| 3 | 09-20-2017 02:37 PM | |
| 1 | 09-20-2017 02:21 PM | |
| 1 | 03-09-2018 03:25 PM | |
| 1 | 03-12-2015 02:06 PM |
| Online Status |
Offline
|
| Date Last Visited |
11-11-2020
02:23 AM
|