|
POST
|
I couldn't resist hacking up a script to see what happened. There were of course a few issues. The MakeFeatureLayer tool would not take a geometry list, so I resorted to using in_memory featureclasses The WITHIN_CLEMENTINI query worked fine interactively in ArcMap but chokes on a polygon with more than a dozen vertices in a script. So a new workaround is needed for that. I test each line as it is generated, and then I could also test its length cumulatively so would only have to test for contains if it is a replacement candidate and I only need to keep the longest current line. This reduced the time per polygon to about 2 seconds. Fairly slow, but it will handle very complex polygons now. Second Hack: #-------------------------------------------------------------------------------
# Name: LongestLine.py
# Purpose: find the longest inside line across a polygon
# Limitations: only single part polygons, no donuts
# Author: kimo
#
# Created: 29/04/2014
# 30/04/2014 changed crossing test to each features
# exclude triangles (use trigonometry for this degenerate case)
# longest line only tested between vertices may need to densify first
#
# Copyright: (c) kimo 2014
# Licence: Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
import arcpy
import sys
import itertools
import datetime
gdb = sys.argv[1]
fcPoly = 'forest'
fcOut = 'fan'
arcpy.env.workspace = gdb
arcpy.env.overwriteOutput = True
begintime = datetime.datetime.now()
desc = arcpy.Describe(fcPoly)
sr = desc.spatialReference
arcpy.env.outputCoordinateSystem = sr
if arcpy.Exists(fcOut):
arcpy.management.Delete(fcOut)
arcpy.management.CreateFeatureclass(gdb,fcOut,"POLYLINE",spatial_reference=sr)
arcpy.management.AddField(fcOut,"ORIGID","LONG")
icur = arcpy.da.InsertCursor(fcOut,['ORIGID','SHAPE@'])
with arcpy.da.SearchCursor(fcPoly,["OBJECTID","SHAPE@"]) as cur:
for row in cur:
gPoly = row[1]
id = row[0]
print "Begin id {} parts {} points {} type {}".format(id,gPoly.partCount,gPoly.pointCount,gPoly.type)
if gPoly.partCount == 1 and gPoly.pointCount > 4: # exclude triangles
starttime = datetime.datetime.now()
lstV = []
for part in gPoly:
for v in part:
lstV.append(v)
oCount = 0
pCount = 0
rCount = 0
maxlength = -1
# use itertools to get unique combinations of pairs
# leave in lines coincident with boundary segments
# because they will always be less than max length if > triangle
for vPair in itertools.combinations(lstV[:-1],2): # exclude last duplicated poly vertex
aLine = arcpy.Polyline(arcpy.Array(vPair),sr)
# dont even bother doing spatial test if already shorter for speed
if aLine.length > maxlength :
if gPoly.contains(aLine): # not lines crossing edges
mLine = aLine
maxlength = aLine.length
pCount+=1
else :
oCount+=1
else:
rCount+=1
print "Max swapped {} Overlap rejects {} Already shorter {}".format(pCount,oCount,rCount)
rec = (id,mLine)
icur.insertRow(rec)
msg = "{} {} time {}".format(id,maxlength,datetime.datetime.now() - starttime)
print msg
arcpy.AddMessage(msg)
del icur
print "done",datetime.datetime.now() - begintime
I added the original Objectid to be able to relate back to each polygon, geometries do not allow me to keep the original ID so I keep it as a python variable.
... View more
04-29-2014
03:25 PM
|
0
|
0
|
11042
|
|
POST
|
That is a ridiculously low performance, especially on a Desktop. You should expect rates in the millions per hour for a single locator. Make sure your reference data is correct and 'perfect'. That is no null or blank values, unique keys and all required fields populated. On ArcServer the peformance will be much slower because you need to send the addresses across a network, process and return the results. In that case performance for bulk geocoding will probably be too slow to be useful. So stick to a network based geocoder for bulk geocoding.
... View more
04-27-2014
06:48 PM
|
0
|
0
|
741
|
|
POST
|
Just generate all the lines by connecting vertex pairs of the polygon and find the longest. Try to avoid iterating around each point combination, think of a process that handles all lines at once. You can leave out the first and last pair of each point iteration because they will be sides, and you only have to do half because the direction does not matter. You can use python collection permutations to get a unique list of pairs. Then turn them into a list of polylines. Because a polygon may not be regular, say a banana shape, then you will need a filter to exclude potential lines crossing the boundary for my simple first stab, but there is a test [ a_poly.overlaps(a_line) ] for that in shape objects. Even better you could test them all in one step, using SelectLayerByLocation(..."WITHIN_CLEMENTINI"...) to select remaining candidates to find the max(length). All the tools will take a geometry list instead of a featureclass so it is fast, in memory and efficient enough to iterate over a large number of polygons.
... View more
04-27-2014
03:37 PM
|
0
|
0
|
11042
|
|
POST
|
How many copies do you expect to sell? I doubt if any password protection is worth the trouble for small obscure scripts such as Model Builder. You can describe the terms in the licence agreement and rely on that to get commercial users to pay. Since you won't get rich by selling scripts, why not just contribute it to the user community as open source and charge for support? Commercial users are likely to be much more interested in your expertise to customise the script for their use than hack it to understand how it works. Consider it advertising for your skills. I have never attempted to hide my scripts, in fact I always provide all the source code for all my work and am still busy. There is always something new to develop even though I still use some scripts that are 20 years old, and I assume my customers may be doing the same. They probably wouldn't be able to do that if they were password protected because no one would remember it.
... View more
04-27-2014
03:06 PM
|
0
|
0
|
550
|
|
POST
|
Maybe we get lulled into placing the output in the same database when the message "All sources must be in the same database" pops up? It does make a tidier workspace to put it with the sources. The locator is 3 files not just one, but in ArcCatalog you only see one entry. But I am now straightened up to place them in a plain old folder. Also note in the fine print that you need the locator in a folder to use multiple threads. Since we all have multicore processors why would we not use them all to get better performance, especially on large datasets?
... View more
04-27-2014
02:55 PM
|
0
|
0
|
3402
|
|
POST
|
Maybe you can make up a legend and convert it into graphic elements? Then ungroup parts so that you can change the legend text. Then you can edit the graphic elements in a Python script to match the parcels. In effect you can fairly easily create your own custom legend if you are using data driven pages. There are several tricks to using graphic elements, give them all a name so you can reference them. Place a lot of copies off the page for each style and just move the required one on to the page as required, say with different numbers of columns. You need to be organised with a list of element names that are easy to iterate over from a list of parcels in the view.
... View more
03-12-2014
02:16 AM
|
0
|
0
|
892
|
|
POST
|
I have just done a couple of tests. Still does it in 10.2, but fixed in 10.2.1 It is all Microsoft's fault for using the escape character as a path separator because the backslash was on the DOS keyboards near the Return (now Enter) key. The rest of the computer world has always used the backslash as an escape for ASCII control characters and something else for path separators. 😮 I could not find a workaround within the dynamic text, but you could run a python script to update the text. Give the text element a name and then use the arcpy.mapping module to access the text and fix it properly. Here is a sample run from the Python window inside ArcMap. Note the element has been named "MXDPath" in the text element properties first. Also since I am using CURRENT, the project must be saved first, Untitled mxds would have a blank filePath. mxd = arcpy.mapping.MapDocument("CURRENT")
mxd.filePath
# u'C:\\GIS\\newFolder\\bugtest102.mxd'
el = arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT","MXDPath")[0]
el.text
# u'Document Path: <dyn type="document" property="path" />'
el.text = mxd.filePath.replace("\\","/") # use the universal path separator supported by Windows for POSIX compliance
el.text
# u'C:/GIS/newFolder/bugtest102.mxd'
mxd.save()
del mxd
You could make this an add-in and trigger it to run when the project is resaved. Note that I am using two backslashes in my replace function because to write a literal backslash you have to escape it with a backslash just like the \n means a newline and \t a tab etc.
... View more
03-12-2014
01:42 AM
|
0
|
1
|
1719
|
|
POST
|
Your SQL query has the quotes mixed up. SQL uses single quotes only for strings, and there is a general kludge to surround invalid field names with double quotes. If the field name has no spaces or punctuation, speclal characters, does not start with a number or an other indecencies in a field name then you can forget the double quotes. This is not the official advice of course, and if you do not have control of field names you may have to add double quotes. SQLexpression = "FIRECODE in ('ARC', 'CEA', 'DAB', 'MGD')" cur = arcpy.SearchCursor(lyr,sQLexpression) Similarly if you really want to have double quotes in the SQL query (why?) instead of escaping them with clumsy escape characters (\) just surround the whole expression with triple quotes of either type. Spaces added for clarity on the triples. SQLexpression = ''' "FIRECODE" in ('ARC', 'CEA', 'DAB', 'MGD') ''' cur = arcpy.SearchCursor(lyr,sQLexpression)
... View more
03-12-2014
12:32 AM
|
0
|
0
|
2336
|
|
POST
|
There are two things you must do: Open the current document in ArcMap using the magic constant 'CURRENT' instead of opening the mxd file. mxd = arcpy.mapping.MapDocument("CURRENT") You also need to run the script in foreground which is a tool property. It is possible to save the CURRENT mxd to a file after modifying it.
... View more
03-11-2014
01:27 AM
|
0
|
0
|
1271
|
|
POST
|
I have never got the locator rebuild tool to work. Ever. I wonder what the tool is for! My solution is to save the result tool that creates the locator in the first place and re-run it again from scratch. I put the results into a small script with a copy-as-a-python-snippet and then use that. It only takes a few minutes to build a new one so you are not saving much time. I write the locator into a local file geodatabase, but this does not package properly for someone else, so I copy the locator into a file before packaging it to copy it into ArcServer.
... View more
03-10-2014
11:47 PM
|
0
|
0
|
4434
|
|
POST
|
A memory leak with only 3000 polygons? I have never experienced that. I cannot think of any process that would not hold all the data in memory and process it in a few seconds. You have said that the script is simple, but your description has a bad feel mentioning a cursor and a tool in the same sentence. If you are running a tool inside a cursor for each feature then I am not surprised you are running out of memory because the tools are not designed to be used that way. Think of the tools like an SQL query - operate on all the features in one operation. You do not loop around a table and run a spatial operation on the featureclass for each record. You should be able to redesign your process to run a tool once. Put it up and let us try to refactor it.
... View more
03-10-2014
11:38 PM
|
0
|
0
|
3304
|
|
POST
|
If you have found an arcobjects solution, it is possible to call arcobjects from Python. Mark Cederholm has had several presentations at Esri conferences on how to do this using comtypes. Maybe this is a last resort?
... View more
02-17-2014
01:25 PM
|
0
|
0
|
1535
|
|
POST
|
Sorry the NTv2http://en.wikipedia.org/wiki/NTv2 is National Transform version 2 developed by the Canadian mapping agency to use a grid to adjust for nonlinear distortions that a datum transform requires. It enables efficient coordinate adjustment that is fast and reliable. It has been adopted by many national agencies to upgrade coordinate systems. It was not relevant to your shift problem because you are not changing datums with a simple shift. I described the latest way to handle projection definitions in ArcGIS that was confusing. At 10.x Esri imbedded all the text file definitions (*.prj) files into a binary dictionary. I have forgotten the magic incantation for the earlier version to create a custom projection. At the worst all you need to do is copy a .prj file and edit with a text editor that handles unix line endings (eg wordpad not notepad). The basic idea still holds, create a custom projection with the offset you require as a false origin.
... View more
02-13-2014
10:32 AM
|
0
|
0
|
2800
|
|
POST
|
Am I missing the point? There was no diagram so maybe I am. Why not dissolve on the selected features before buffering? Then the adjacent boundaries will not be buffered. Are the adjacent polygons not touching? In that case there are other tools that can jump the gap in the Cartographic Toolset. eg AggregatePolygons Tool
... View more
02-12-2014
12:52 AM
|
0
|
0
|
7400
|
|
POST
|
Use a custom projection with a false origin of 360 degrees. We have just been doing shifts in NZ after major earthquakes which have all features moving different amounts over a lot of featureclasses. Too many for rubber sheeting to handle. The solution was to create a transform grid using the NTv2 tools and run a pseudo-projection by defining a custom projection. After the projection using transform we redefine the projection back to the original. In your case it is easier because you don't need a transform, just a shift. Get out a copy of the projection (from properties of a featureclass in ArcCatalog) and hit the copy_and_modify button. Add a false origin of 360 degrees for your custom projection. Then reproject to the new projection and afterwards override the new projection by redefining it back to the original projection. Merge the two.
... View more
02-12-2014
12:34 AM
|
0
|
0
|
2800
|
| Title | Kudos | Posted |
|---|---|---|
| 1 | 2 weeks ago | |
| 1 | 03-11-2023 03:54 PM | |
| 1 | 09-15-2024 10:32 PM | |
| 1 | 03-12-2026 01:10 AM | |
| 1 | 03-13-2026 08:30 PM |
| Online Status |
Offline
|
| Date Last Visited |
2 weeks ago
|