|
POST
|
I have found the interactive Merge tool with a fieldmap is a good substitute for a join. It works very efficiently in ArcMap but not in a script. I just reopen the results in the saved MXD and rerun the tool again. It is so much faster that its worth the inconvenience once a month on a pair of files containing 2M records. Maybe one day it will work in a script? I did solve another problem with multiple joins never finishing by partitioning the task into 12 parts because there happened to be a key I could use for subsets. Each part ran in less than 10 minutes for under an hour of processing. If you have a way of selecting subsets that may avoid exceeding your memory limits and still use the python dictionary trick. I limit my count to about a million records and it seems to be OK.
... View more
11-08-2012
01:05 AM
|
0
|
0
|
5247
|
|
POST
|
Fails for me every time in a script, but not interactively. A complete pain. Even if I run it successfully interactively, copy-as-python-snippet and repeat in a script it fails to populate the joined fields. There is no error message, just blank (not null) fields. Running 10.0 SP5 but a test on 10.1 SP1 had the same result. My (terrible) workaround is to create a dictionary of the joined feature content, add the fields and then run an updatecursor with the dictionary. Alternatively I just open the results the next month and run it again while I surf the net watching the election.....
... View more
11-08-2012
12:49 AM
|
0
|
0
|
2610
|
|
POST
|
You don't have to add fields one at a time. Use the FeatureClasstoFeatureClass tool and a Fieldmap. This does not require programming and can be included in your model. You can define new fields there. While you are there, maybe you can merge the source data so you would not have to calculate the fields either. It is very fast if all you were doing was doing a join and calculate which should be avoided because it is very slow and often leaves the fields blank.
... View more
11-08-2012
12:40 AM
|
0
|
0
|
2181
|
|
POST
|
I have some that take hours or even days That would seem unreasonable in my experience. No single geoprocessing step or model should take days. "If a process takes longer than a cup of coffee, then interrupt the process and find a better way." Is my rule. Things to check: Do not run processes on a desktop with data across a network Ensure that all key fields are indexed Use filegeodatabases, not shapefiles, but if you must, index the shapefiles If you are running geoprocessing tools inside a loop of features, recode it so it is a single pass like using SQL set commands. Optimise your disk, have 30% free space, plenty of virtual memory. Make your scratch workspace a filegeodatabase located on a local disk, this is critical for modelbuilder Abandon modelbuilder and recode the processes in Python scripts Avoid known slow tools, such as calculating field values with joined tables Use selections to minimise the unnecessary movement of data Partition the process to keep within memory limits Split the task and run parallel tasks on a multiprocessor machine Add a faster disk with a RAID0 array Maybe look for some different software to do the task, eg FME, numpy, PIL, gdal.
... View more
11-08-2012
12:28 AM
|
0
|
0
|
3456
|
|
POST
|
You have not specified how you are reading in the text file. Have you had a look at the schema.ini specification? This is a Microsoft standard for defining the schema for reading in text files. It works really well for loading text files into a geodatabase table because you can rename, set input widths, types to make the simple table to table tool work very well. This avoids a lot of complex coding opening a file and unpicking the components. If you have a table in a database try exporting it as a txt file and at 10.x a schema.ini file will be generated for you. Then keep the file with the source, edit the name and hey presto a schema that works. Here is an example where I actually write the schema after dumping out the table as a text file. This is much faster than an InsertCursor at 10.0. Maybe the new da module makes it faster, but it is still a very easy way of defining a schema. ..... f2 = open("e:/crs/currentgdb/schema.ini","w") f2.write("""[owners.csv] Format=CSVDelimited ColNameHeader=True col1=TTL_TITLE_NO Text Width 20 col2=OWNERS Text Width 120""") f2.close() print "uses schema.ini" arcpy.conversion.TableToTable(csvOwners,ws,"owners")
... View more
11-07-2012
11:34 PM
|
0
|
0
|
2136
|
|
POST
|
Your description appears to be a bit loose. Do you mean you want to calculate the bearing of the segments between bends, or the final heading of the end of the pipeline, or the included angle? Maybe you want to add a bend in the line given the offset implied by the joint name? The math module has a function atan2() that returns a 360 degree bearing given the delta_x and delta_y. There are also functions to convert radians to degrees and all the trig functions.
... View more
10-13-2012
12:54 AM
|
0
|
0
|
1011
|
|
POST
|
I have a metre of bookshelf of Python books, but there is one that is really dog-eared. It is out of print, but can be picked up on Amazon. Be careful not to get the second edition, it is a completely different book by a different author that is not the same format and is not useful. Some reprints are very poor quality too. "Python" by Chris Fehily, Peachpit Press ISBN 0-201-74884-3. Most books I have are far too technical for what is needed for geoprocessing. http://www.amazon.com/Python-Chris-Fehily/dp/0201748843/ref=cm_cr_pr_pb_t The simple layout, examples and explanations are just the right level for writing scripts. True, it was written for Python 2.2 but there are other books for version 2.7 or 3000. I bought "The Quick Python Book", Ceder, Manning which covers Python 3 and have not found it useful.
... View more
10-13-2012
12:42 AM
|
0
|
0
|
493
|
|
POST
|
I thought that this was fixed in 10.1 to use the standard install. I prefer python to be in C:\Python27 so that all other modules have this expected location. Some things such as the python help seem to be hardcoded for this location, not c:\Python27\ArcGIS10.1 All you have to do is save the Desktop10.1.pth file in site-packages, and the other modules numpy and matplotlib, uninstall python, delete the pesky directory and reinstall python2.7.x from the standard distribution. Then copy back the modules and the pth file which locates the arcpy modules. Each time you install a service pack, rerun the python install to repair the paths back to normal. I have never had any problem with any dot release of the main version eg 2.6.x or 2.7.x which should only be bug fixes. You may not get the bug fixed if you run the script in-process because the python version is imbedded. While you are there, install the Pythonwin package which gives you extra modules to win32 API, COM interfaces to MS Office applications, message boxes and even the registry. A new location? Be careful what you do about that. You cannot have any spaces in the path. Maybe you are an 'administrator' who thinks that c:\Program Files (x86) would be tidier? You will condemn everyone to inexplicable crashes and error messages.
... View more
08-08-2012
05:09 PM
|
0
|
0
|
961
|
|
POST
|
If you are in a Python script already I recommend that you do not use CalculateField, use a cursor. CalculateField is useful in ModelBuilder and interactive table manipulation but not in a script. It is no faster than a cursor (especially at 10.1 da module) because CalculateField simply wraps a cursor around the expression. With a cursor you can: Trap errors with a try/except block Test for unexpected data and handle it properly (eg null values) Debug easily by adding print statements Write the expressions without complicated syntax (no !wrappers! or [access_speclal] brackets) Do multiple field calculations at once Use a dictionary to update values to avoid a join
... View more
08-08-2012
04:49 PM
|
1
|
0
|
1090
|
|
POST
|
Yes, don't bother with the string. Create a SpatialReference object and set the Factory Code or the Name.
sr = arcpy.SpatialReference()
sr.factoryCode = 32759 # case sensitive
sr.create() # essential step that is non-intuitive!
print sr.name # to check
At 10.1 the library of strings for standard projections has disappeared! They are now built in to ArcGIS. The factorycode is the most reliable but hardest to find, look up the projected_coordinate_systems.PDF in the root ArcGIS documentation directory. I can never remember what name property is required to be set to make it work.
... View more
08-08-2012
04:39 PM
|
0
|
0
|
7021
|
|
POST
|
Another example is with time zones. I like pytz and I'm used to using it, but ArcGIS ships with its own time zone database and its own way of doing things with the time-aware stuff that came with 10.0, so instead of bundling pytz with its own parallel time zone database I went ahead and integrated the ArcGIS way of doing things as arcpy.time, which had the additional benefit of the smart time delta class that uses ArcGIS' time APIs. This then also negated any specific need for dateutil, which is another library I really like. But arcpy.time handles standard Python datetime and tzinfo objects so it all interoperates just fine with everything else, and the skills you learn either from arcpy or pytz/dateutil transfer without a hitch. When I backported GPX to Features in 10.1 I found that the time conversion tools do not work even at SP5 so I have to go back to dateutil.
... View more
08-08-2012
02:18 PM
|
0
|
0
|
2824
|
|
POST
|
Shipping on the CD is not very helpful. Do you really mean that these modules should be installed with Python? I have a current problem writing some tools to generate Excel spreadsheets as final reports. Easy to do using the win32com module that comes with win32all or Pythonwin that is included with the standard install but not installed with ArcGIS and Python. The result is that scripts will work in a development environment but cannot be deployed because the modules are not organisation wide. Is every IT administration department as difficult to deal with as my students encounter? They are very unhelpful installing service patches, python modules, customising Window to include file extensions, and dozens of other settings that impact on ArcGIS users. I suggest various system default settings to optimise ArcGIS and they just roll their eyes. A very serious performance hit is placing the default.gdb in the user's profile area which is not local, but across the network. So rather than worry about some optional and less used modules, I would rather see the defaults loaded with the software to tackle intransigent and unskilled administrators. I suppose they can always reverse them to 'comply with Policy'. As to having the lowest common denominator in Linux not having win32 capability? That is not an argument for not having it for the majority that do use Windows applications. After all ArcMap only works on Windows, by that argument we should not have that because it does not run on Linux.
... View more
08-08-2012
02:04 PM
|
0
|
0
|
2824
|
|
POST
|
Loading dates is always tricky. Do you have a schema.ini file to define the schema of your csv file? That will help a lot. You will also have to have your date defined in a format that Microsoft can understand. You specify the DateTimeFormat in the schema.ini file. If that is not working, then you will need to run a script that reads in the date as a string field, parse it to a python datetime object and then write it out to a date field. I hope you are not using shapefiles, because that can only handle a date without time. See the GPX to Featureclass conversion tool in 10.1 for an example. Unfortunately the tool to do this conversion does not work in 10.1, so I wrote a Python snippet: from datetime import datetime from dateutil import tz ## xxxxxxxxxxxxxx snip xxxxxxxxxxxxxxxxx cur = arcpy.InsertCursor(outFC,srWGS84) # Loop over each point in the tree and put the information inside a new row # for index, trkPoint in enumerate(GeneratePointFromXML(tree)): if trkPoint.asPoint() is not None: ## rowsDA.insertRow([trkPoint.name, trkPoint.desc, trkPoint.gpxtype, trkPoint.t, ## float(trkPoint.z), float(trkPoint.x), float(trkPoint.y), float(trkPoint.z)]) row = cur.newRow() row.Name = trkPoint.name row.Descript = trkPoint.desc row.Type = trkPoint.gpxtype row.DateTimeS = trkPoint.t # convert UTC to nz using Python because ArcGIS tool crashes at 10.0 utc = datetime.strptime(trkPoint.t,'%Y-%m-%dT%H:%M:%SZ').replace(tzinfo=from_zone) # Tell the datetime object that it's in UTC time zone since # datetime objects are 'naive' by default # Convert time zone back for ArcGIS local = utc.astimezone(to_zone).replace(tzinfo=None) row.DateTime = local row.Elevation = float(trkPoint.z) row.shape = arcpy.Point(float(trkPoint.x), float(trkPoint.y), float(trkPoint.z)) cur.insertRow(row) else: badPt +=1
... View more
08-07-2012
03:39 PM
|
0
|
0
|
640
|
|
POST
|
If you were a keen scripter you could probably write a much faster change detector in Python using the SET datatype. You have to load the sets of course, using a SearchCursor and creating a dictionary using the key. This takes a few seconds, but the result is definitely in memory. The set comparison takes milliseconds. Then you have to painfully write out the difference records using a layer. Depending on how many changes this can take a few minutes, or in your case, an update script can loop through the list. How large are your tables, and how many differences do you have? For an idea, I do a comparison of two featureclasses of polygons containing 2.45 million records. Last month there were 11616 changes, 1516 deletions and 2270 new records which are written out to new featureclasses. It takes around 40 minutes to run. Here is another example. When creating layer definition views you MUST index by the key for the SQL query to work in human time. #NAME: buildDeltaLYRs.py
#AUTHOR: Kevin Bell
#EMAIL: [email protected]
#DATE: 20071207
# edited by KimO
#PURPOSE: create adds/deletes layer files by comparing 2 point
# feature classes shape and attributes. If the shape has
# not changed, but any of the attributes have, the feature
# will show as a delete, and an add.
#RECOMMENDED SYMBOL: adds.lyr = green plus, deletes.lyr = red X
# (this allows for nice stacking)
#NOTE: __buildDict method has hard coded primary key and attribute names.
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
def buildDeltaLayers(inFC1, inFC2):
'''build an adds and deletes lyr for a given chrono FC '''
d1 = __buildDict(inFC1)
d2 = __buildDict(inFC2)
compareList = __valuesChanged(d1,d2)
__makeLYR(inFC1, compareList, "deletes")
print "...deletes.lyr is complete."
d1 = __buildDict(inFC2)
d2 = __buildDict(inFC1)
compareList = __valuesChanged(d1,d2)
__makeLYR(inFC2, compareList, "adds")
print "...adds.lyr is complete."
def __valuesChanged(dict1, dict2):
'''get a list of keys from one dict if a cooresponding dict's values are different '''
outList = [key for key in set(dict1.keys() + dict2.keys()) if dict1.get(key) != dict2.get(key)]
return outList
def __buildDict(inputFC): #-----BEWARE OF HARDCODED PRIMARY KEY AND ATTRIBUTES BELOW!!!!!
'''Build a dictionary of the primary key, and it's fields'''
d = {}
cur = gp.SearchCursor(inputFC)
row = cur.Next()
while row:
d[row.FP] = [row.Shape.Centroid, row.LUMENS, row.WATTS, row.TYPE, row.OWNER, row.RATE]
row = cur.Next()
del cur
return d
def __makeLYR(fc, inList, outLyrName):# BEWARE OF HARDCODED PRIMARY KEY BELOW
'''given a list, return a LYR file'''
wc = str(tuple(inList))
whereclause = "FP IN " + wc # <----IF DATA ISN'T FILE GDB, YOU MAY NEED QUOTES/BRACKETS
gp.MakeFeatureLayer_management (fc, outLyrName, whereclause)
gp.SaveToLayerFile_management (outLyrName, outLyrName +".lyr")
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
print "----------- building adds/deletes layer files ----------------"
import arcgisscripting, time
startTime = time.clock()
gp = arcgisscripting.create()
gp.workspace = r"F:\Temp\streetLightCompare.gdb"# <----BEWARE OF HARDCODED WORKSPACE !!!!!
gp.OverwriteOutput = 1
#o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o
buildDeltaLayers("SEP07","OCT07") #<-------------------BEWARE OF HARDCODED FEATURE CLASSES !!!!!
#o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o
print "Your adds/deletes lyr's are in:"
print str(gp.workspace)
del gp
print "--------------------------------------------------------------------------------------------"
stopTime = time.clock()
elapsedTime = stopTime - startTime
print "elapsed time = " + str(round(elapsedTime, 1)) + " seconds"
... View more
08-07-2012
03:28 PM
|
0
|
0
|
1430
|
|
POST
|
It has taken me years to work out how to use fieldmapping(s) objects, so 4 hours is just a start. Here is a snippet of code to help. What I am doing is merging a library of coverage tiles (remember them? still useful) My most reliable way to get it to work is to create a fieldmappings object from the source and then delete fields that I don't want from a list. You could reverse the logic and delete all fields except the ones you do want. Note that when I merge I have to add each input source fieldmap, just having the output names once will not do! The alternative is to use the CopyFeatures tool interactively and copy the result as a python snippet for a static version. It produces a string equivalent of a fieldmap which is not ideal, and really hard to edit. I do this when scripting fails and I have to run the tool interactively. I just get out last month's result, open it as a dialog and rerun. start = datetime.datetime.now() droplst = ["xxx","rings_ok","rings_nok","lpoly#","rpoly#", "fnode#","tnode#","length","area","perimeter", "multi","case1","case2","case3", "$SCALE","$ANGLE","PAR-ID","OTHR#","OTHR-ID",] sr = arcpy.SpatialReference() sr.factoryCode = 2193 # NZTM sr.create() arcpy.env.outputCoordinateSystem = sr ws = "e:/lib/nztm/tile" arcpy.env.workspace = ws+"/t1001/data" dictCov = { "add_all":"address/point", "emf":"emf/point", "fename":"fename/point", "hydro":"hydro/polygon", "titledif":"titledif/point", "up":"up/point", "uparcel":"uparcel/point" } l # lstFC = ["uparcel"] for outFCName in lstFC : cov = dictCov[outFCName] fds = cov.split("/")[0] print "Processing",outFCName,cov arcpy.AddMessage(outFCName) if not arcpy.Exists(gdb+"/"+outFCName) : lstSrc = [] for ld in range(12) : tile = "t%d" % (ld + 1001) srcCov = ws+"/"+tile+"/data/"+cov if arcpy.Exists(srcCov) : lstSrc.append(srcCov) else : arcpy.AddError(srcCov+'not found') print srcCov,"NOT FOUND" # Create FieldMappings object to manage merge output fields fieldMappings = arcpy.CreateObject("FieldMappings") # Add all fields from sources for src in lstSrc: fieldMappings.addTable(src) # hide fields print " removing ", for fld in [fds+"#",fds+"-ID"] + droplst: #include cov# and cov-id pos = fieldMappings.findFieldMapIndex(fld.upper()) # print pos,fld if pos >= 0 : fieldMappings.removeFieldMap(pos) print fld.upper(), print print " keeping ", for x in range(fieldMappings.fieldCount): print fieldMappings.getFieldMap(x).outputField.name, print arcpy.Merge_management(lstSrc, gdb+"/"+outFCName,fieldMappings ) print outFCName,"merged" arcpy.AddMessage(outFCName+" merged") else : print "Skipping load of",outFCName arcpy.AddWarning("Skipping load of "+outFCName)
... View more
08-07-2012
03:00 PM
|
0
|
0
|
1364
|
| Title | Kudos | Posted |
|---|---|---|
| 3 | yesterday | |
| 1 | 3 weeks ago | |
| 1 | 03-11-2023 03:54 PM | |
| 1 | 09-15-2024 10:32 PM | |
| 1 | 03-12-2026 01:10 AM |
| Online Status |
Offline
|
| Date Last Visited |
yesterday
|