|
POST
|
Sounds as if your subprocess script is failing immediately after it is called. It might be useful for it to write some sort of log file so you can see how far it gets. for calling scripts from a "parent" script I have always used the old fashioned os.spawnv. For example: rootDir = r"D:\temp"
script = r"\\snarf\div_lm\script.py"
dateTimeStamp = time.strftime('%Y%m%d%H%M%S')
for path in sys.path:
if os.path.exists(os.path.join(path, "python.exe")) == True:
pythonExePath = os.path.join(path, "python.exe")
parameterList = [pythonExePath, script, rootDir, dateTimeStamp[:8]]
os.spawnv(os.P_WAIT, pythonExePath, parameterList)
#or to not wait for the spawned process to finish...
os.spawnv(os.P_NOWAIT, pythonExePath, parameterList) If you figure out how to get the subprocess module to use a "no wait" method, can you post it? I never found a way to do it...
... View more
07-14-2010
08:29 AM
|
0
|
0
|
1828
|
|
POST
|
Locks often remain for the duration of an ArcMap or ArcCatalog session, even though you aren't currently looking/accessing the FGDB. Also, the gp object does the same thing! Make sure you don't have any ArcWhatever or python.exe running (that are looking or ever looked at one time at the FGDB), and then run your script.
... View more
07-08-2010
03:25 PM
|
0
|
0
|
2377
|
|
POST
|
The fishnet tool is basically a way to make a polygon layer that looks like, well, a fishnet that has adjacent rectangualr (or square) polygons (like a checkerboard). I do a lot of large-dataset processing in Python, and the Fishnet tool is always my step #1 in being able to break up large datasets to get them to work with the ESRI tools. For example: # Description
# -----------
# This creates an area of interest feature class, a few buffers of the area of interest, and a fishnet polygon to use for tile-by-tile psudo-parallel processing.
# Written for Python version: 2.5.1 (PythonWin)
# Written for ArcGIS version: 9.3.1 (should be forward compatible).
# Author: Chris Snyder, WA Department of Natural Resources, chris.snyder(at)wadnr.gov
# Date created:
# Last Updated: 20090923, csny490
try:
#Import system modules and creates the Geoprocessor object (the v9.2 way)
import sys, string, os, glob, shutil, time, traceback, arcgisscripting
gp = arcgisscripting.create(9.3)
#Defines some functions
def showGpMessage():
try:
print >> open(logFile, 'a'), "\n" + gp.GetMessages() + "\n"
print "\n" + gp.GetMessages() + "\n"
except:
pass
def showPyMessage():
try:
print >> open(logFile, 'a'), str(time.ctime()) + " - " + str(message)
print str(time.ctime()) + " - " + str(message)
except:
pass
#Sets a date/time stamp variable in the format yyyymmddhhmmss (example: 20060330132345)
dateTimeStamp = time.strftime('%Y%m%d%H%M%S')
#Specifies the root variable, makes the logFile variable, and does some error checking...
root = sys.argv[1]
if os.path.exists(root)== False:
print "Specified root directory: " + root + " does not exist... Bailing out!"
sys.exit()
scriptName = sys.argv[0].split("\\")[len(sys.argv[0].split("\\")) - 1][0:-3] #Gets the name of the script without the .py extension
logFile = root + "\\log_files\\" + scriptName + "_" + dateTimeStamp[:8] + ".log" #Creates the logFile variable
if os.path.exists(root + "\\log_files") == False:
os.mkdir(root + "\\log_files")
if os.path.exists(logFile)== True:
os.remove(logFile)
message = "Deleting log file with the same name and datestamp... Recreating " + logFile; showPyMessage()
workspaceDir = "overlay_index"
if os.path.exists(root + "\\" + workspaceDir)== False:
message = "Creating layer directory: " + root + "\\" + workspaceDir; showPyMessage()
os.mkdir(root + "\\" + workspaceDir)
#Attempts to check out the highest grade license available...
if gp.CheckProduct("ArcInfo") == "Available":
gp.SetProduct("ArcInfo")
elif gp.CheckProduct("ArcEditor") == "Available":
gp.SetProduct("ArcEditor")
elif gp.CheckProduct("ArcView") == "Available":
gp.SetProduct("ArcView")
else:
messsage = "No ArcGis licesnses are available... Exiting script!"; showPyMessage(); sys.exit()
message = "Selected an " + gp.ProductInfo() + " license"; showPyMessage()
#Sets some environment variables
gp.overwriteoutput = True
gp.loghistory = False
sr = 'PROJCS["NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",45.83333333333334],PARAMETER["Standard_Parallel_2",47.33333333333334],PARAMETER["Latitude_Of_Origin",45.33333333333334],UNIT["Foot_US",0.3048006096012192]]'
gp.OutputCoordinateSystem = sr
gp.XYTolerance = "0.001 METERS"
gp.XYResolution = "0.0005 METERS"
gp.workspace = root + "\\" + workspaceDir
#Process: Creates a FGDB to hold all the spatial layers
fgdbName = "overlay_index"
fgdbPath = gp.workspace + "\\" + fgdbName + ".gdb"
gp.CreateFileGDB_management(gp.workspace, fgdbName); showGpMessage()
#Process: Makes a FC of the PLS polys we want
plsFC = root + "\\gis_layers\\pls.gdb\\pls"
plsFL = "pls_feature_layer"
gp.MakeFeatureLayer_management(plsFC, plsFL, "SUR_OWN_CD > 0 OR TIM_OWN_CD > 0", "", ""); showGpMessage()
dnrMgmtAllFC = fgdbPath + "\\dnr_mgmt_all"
gp.Dissolve_management(plsFL, dnrMgmtAllFC, "", "", "SINGLE_PART"); showGpMessage()
#Process: Creates a fishnet layer (the basis of the index polygons)
fishnetLine1FC = fgdbPath + "\\fishnet_line1"
dsc = gp.describe(dnrMgmtAllFC)
xMin = dsc.extent.xmin
yMin = dsc.extent.ymin
xMax = dsc.extent.xmax
yMax = dsc.extent.ymax
bufferWidth = 1 #Buffer width expands the the extent of the fishnet by that many map units (to make sure it extends past DNR lands a bit)
numberOfRows = 6
numberOfColumns = 6
gp.CreateFishnet_management(fishnetLine1FC, str(xMin - bufferWidth) + " " + str(yMin - bufferWidth), str(xMin - bufferWidth) + " " + str(yMin), "0", "0", numberOfRows, numberOfColumns, str(xMax + bufferWidth) + " " + str(yMax + bufferWidth), "NO_LABELS", ""); showGpMessage()
fishnet1FC = fgdbPath + "\\fishnet_poly1"
gp.FeatureToPolygon_management(fishnetLine1FC, fishnet1FC, "", "", ""); showGpMessage()
gp.Delete_management(fishnetLine1FC, ""); showGpMessage()
#Process: Creates a 2nd fishnet layer (that is more dense than the original)
fishnetLine2FC = fgdbPath + "\\fishnet_line2"
densificationFactor = 12 #this densification factor ABSOLUTELY NEEDS to be a multiple of 2 (e.g. 2,4,6,8,etc)
gp.CreateFishnet_management(fishnetLine2FC, str(xMin - bufferWidth) + " " + str(yMin - bufferWidth), str(xMin - bufferWidth) + " " + str(yMin), "0", "0", str(numberOfRows * densificationFactor), str(numberOfColumns * densificationFactor), str(xMax + bufferWidth) + " " + str(yMax + bufferWidth), "NO_LABELS", ""); showGpMessage()
fishnet2FC = fgdbPath + "\\fishnet_poly2"
gp.FeatureToPolygon_management(fishnetLine2FC, fishnet2FC, "", "", ""); showGpMessage()
gp.Delete_management(fishnetLine2FC, ""); showGpMessage()
#Process: Updates fishnet1FC with fishnet2FC for specified areas
hcpUnitFC = root + "\\gis_layers\\hcpunit.gdb\\hcpunit"
oesfFC = fgdbPath + "\\oesf"
gp.Select_analysis(hcpUnitFC, oesfFC, "HCPUNIT_ID = 1", "", ""); showGpMessage()
fishnet1FL = "fishnet_1_feature_layer"
gp.MakeFeatureLayer_management(fishnet1FC, fishnet1FL, "", "", ""); showGpMessage()
gp.SelectLayerByLocation_management(fishnet1FL, "INTERSECT", oesfFC, "", "NEW_SELECTION"); showGpMessage()
fishnet2FL = "fishnet_2_feature_layer"
gp.MakeFeatureLayer_management(fishnet2FC, fishnet2FL, "", "", ""); showGpMessage()
gp.SelectLayerByLocation_management(fishnet2FL, "HAVE_THEIR_CENTER_IN", fishnet1FL, "", "NEW_SELECTION"); showGpMessage()
gp.Delete_management(oesfFC, ""); showGpMessage()
fishnetFC = fgdbPath + "\\fishnet_poly"
gp.Update_analysis(fishnet1FC, fishnet2FL, fishnetFC, "BORDERS", ""); showGpMessage()
#Process: Make a line FC version of fishnetFC
fishnetLineFC = fgdbPath + "\\fishnet_lines"
gp.FeatureToLine_management(fishnetFC, fishnetLineFC, "", "NO_ATTRIBUTES"); showGpMessage()
#Process: Unions everything together
union1FC = fgdbPath + "\\union1"
gp.Union_analysis(fishnetFC + ";" + dnrMgmtAllFC, union1FC, "ONLY_FID", "", "GAPS"); showGpMessage()
#BLAH BLAH BLAH - I had to trim this part out to post it to the forum (stupid 10,000 charater limit!!!)
#Process: Compacts the PGDB
gp.Compact_management(fgdbPath); showGpMessage()
#Process: Creates a .txt file that indicates to other scripts how many tiles there are (this is used so that the run_union_master script doesn't have to use the gp and waste memory)
print >> open(root + "\\log_files\\" + "therearethismanytiles_" + str(numberOfTiles) + ".txt", 'a'), "therearethismanytiles_" + str(numberOfTiles)
message = "ALL DONE!"; showPyMessage()
print >> open(root + "\\log_files\\" + scriptName + "_isalldone.txt", 'a'), scriptName + "_isalldone.txt"
except:
print >> open(root + "\\log_files\\" + scriptName + "_bombed.txt", 'a'), scriptName + "_bombed.txt"
message = "\n*** LAST GEOPROCESSOR MESSAGE (may not be source of the error)***"; showPyMessage()
showGpMessage()
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()
message = "\n*** PYTHON LOCAL VARIABLE LIST ***"; showPyMessage()
variableCounter = 0
while variableCounter < len(locals()):
message = str(list(locals())[variableCounter]) + " = " + str(locals()[list(locals())[variableCounter]]); showPyMessage()
variableCounter = variableCounter + 1
showPyMessage()
... View more
07-08-2010
09:54 AM
|
0
|
0
|
2643
|
|
POST
|
In v9.3.1 ArcInfo level, the Clip_managment tool does indeed have an option for using a vector featureclass as a psudo-mask (see the "in_template_dataset" parameter). This is tool to use when you have a multi-band image, since the output can be any ESRI-supported format (not just GRID like the ExtractByMask tool). Check it out...
... View more
07-08-2010
08:47 AM
|
0
|
0
|
1426
|
|
POST
|
If workstation fails (try the regionbuffer tool), you might try this in ArcGIS: 1. Convert your lake polys to a line (FeatureToLine tool) 2. Buffer the lines 20m (on both sides, round ends). If the shoreline buffer fails try: a. Build a fishnet polygon (CreateFishnet tool and then FeatureToPolygon) of sufficient # of cells to break the large contiguous line features into smaller segments. This might take a couple tries to get right. b Buffer the "chopped up" lines. 3. Erase the buffer from the lake poly. If the erase fails try: a. Intersect the lake polys with the fishnet you created for the buffers b. Try the intersect again with the "chopped up" lake polys Generally, the issue is that ArcGIS chokes when you try to do some sort of geoprocessing operation with a feature that has a kajillion vertices. The reason being that all the feature's x,y, pairs need to get loaded into memory in order to do the operation, and often the number of vertrices in a complex feature exceed the limits of the ArcGIS algorithm. The solution is to break your input features into small enough pieces to allow the algorithms to run without maxing out the memory limitations.
... View more
07-08-2010
08:17 AM
|
0
|
0
|
2643
|
|
POST
|
Try: 1. Breaking the input FC into singlepart geometry using the MultipartToSinglepart tool 2. Running the RepairGeometry tool. Are there any "errors"? 3. Then run the negative buffer. Does it work then? The issue of course is the extreemly complex geometry of the lake polygons. If the buffer tools fails on the original FC, maybe try writing a script that will loop through the features and make a buffer for each feature in a try/except block. That way you will be able to figure out the offending features (maybe it's only a few of them?).
... View more
07-07-2010
09:43 AM
|
0
|
0
|
2643
|
|
POST
|
I know how to write a script, but don't have much experience with doing batch processing with ModelBuilder (although I'm sure it can be done there). A script in ArcGIS v9.3 would look something like this: import arcgisscripting, shutil
gp = arcgisscripting.create(9.3)
gp.workspace = r"C:\temp\this_is_where_my_orig_tif_files_are"
tifList = gp.listrasters("","TIFF")
for tif in tifList:
gp.ExportRasterWorldFile_management (tif)
... View more
07-06-2010
12:28 PM
|
0
|
0
|
603
|
|
POST
|
Can't do it in v9.3, but now in v10.0 there is a zoomToSelectedFeatures() method as part of the arcpy.mapping module. This would probably have to be in script format (not modelbuilder) to work though....
... View more
07-06-2010
11:28 AM
|
0
|
0
|
2885
|
|
POST
|
check out the SpatialReference parameter in the search cursor. Unfortunately you cannot directly convert coordinates using gp.describe, so you would either have to: 1. Using the SpatialReference parameter in the search cursor, loop through all the features keeping track of the xMin, yMin, xMax, and yMax extents of each feature, and then compute the xMin, yMin, xMax, and yMax for the entire layer. 2. Using gp.describe, get the layers extent, construct a polygon feature for it, project it (using the project tool), and then reread the coordinates.
... View more
07-06-2010
09:28 AM
|
0
|
0
|
524
|
|
POST
|
1) $ARCHOME is a environment variable for ArcInfo Workstation. If you are using a Windows machine, you can find what path this variable points to in many ways. But probably the easiest is to: a. Go to your desktop, and right click on My Computer > Properties > Advanced (tab) > Environment Variables. b. Look under the System Variables, and you should see one for ARCHOME (probably something like C:\Program Files\.... 2) Okay, I get it. Question: Did your original .tif files have accompanying .tfw files? If so, all you need to do is copy those into the location where your new photoshopped .tif files are, and make sure the have the same prefix name (for example: photoshop_test4.tfw goes with photoshop_test4.tif). If your originals .tif files (the ones you used as input to your Photoshop process) did not have .tfws (thus indicating they were either geotiffs - or were not projected in the 1st place?!?), it will be a bit harder, but still pretty easy. Take a look at the ArcGIS tool called " Export Raster World File (Data Management)". Basically you would have to write a script that would look through all your original .tif files, run the tool on them, and then write the output .tfw files to the path of your new post-photoshopped images. There are harder ways to do this as well, but let me know if either of these two methods works.
... View more
07-06-2010
08:38 AM
|
1
|
0
|
4646
|
|
POST
|
Seems like your value99 variable is numeric, right? Maybe as simple as converting value99 to a string so its "concatenatable" with the SOMA exp? gp.SingleOutputMapAlgebra_sa("CON( " + raster + " >= " + str(value99) + ", 99, 0)", Output)
... View more
06-30-2010
01:56 PM
|
0
|
0
|
594
|
|
POST
|
Check out the Clip_managment tool (not the Clip_analysis one). Be sure to set the snap raster in the Envr. Settings, otherwise the pixel allignment may differ from the original.
... View more
06-29-2010
09:30 AM
|
0
|
0
|
1426
|
| Title | Kudos | Posted |
|---|---|---|
| 1 | 08-29-2024 08:23 AM | |
| 1 | 08-29-2024 08:21 AM | |
| 1 | 02-13-2012 09:06 AM | |
| 2 | 10-05-2010 07:50 PM | |
| 1 | 02-08-2012 03:09 PM |
| Online Status |
Offline
|
| Date Last Visited |
08-30-2024
12:25 AM
|