Latest Contributions by JoshuaBixby

Overall, I think the ArcPy support for Z and M coordinates is poor, quite poor.  I wish I had a better response, but I just hit walls/bugs/design limitations around every corner. I quickly revisited arcpy.AsShape and arcpy.FromWKT.  Although the Help for FromWKT doesn't state anything one way or another about Z or M coordinates, the function throws an error when you try to use Z or M coordinates.  I have been able to use arcpy.ArcSDESQLExecute to pass Z and M WKT coordinates to SQL Server Express, basically bypass ArcPy and relying on the SQL Server engine to handle creating the shapes correctly.  Even though it has worked for me, and I can view the results in ArcMap, all of the regular and Data Access cursors seem to do nothing with them.  It might just be ArcPy cursors can't cope with Z or M coordinates, need to investigate more. Good luck, let us know if you have any success. UPDATE:  I had to strike through some of my earlier comments.  On the WKT comment, it is still the case that the Help doesn't explicitly mention Z or M coordinates; however, arcpy.FromWKT was throwing an error because I wasn't giving it OGC compliant WKT representations.  I am used to MS SQL Server's WKT implementation, which is slightly different than the OGC specifications, so what I would pass SQL Server was failing with ArcPy. On the ArcPy cursors not working with ZM polygons created through SQL Server and arcpy.ArcSDESQLExecute, I had a goofy projection discrepancy that was causing the problems.  Once I went back and straightened out the projections, I was able to use cursors with the newly created polygons. ... View more

I agree that Python object equivalence/equality likely isn't the best approach for comparing spatial objects primarily because I haven't seen any geospatial Python packages where the various geometry classes' __eq__() methods are spatially aware.  With ArcPy, they are not. >>> pl1 = arcpy.Polyline( ...     arcpy.Array( ...         [arcpy.Point(0, 0), arcpy.Point(1, 1)] ...     ) ... ) ... >>> pl2 = arcpy.Polyline( ...     arcpy.Array( ...         [arcpy.Point(0, 0), arcpy.Point(0.5, 0.5), arcpy.Point(1, 1)] ...     ) ... ) ... >>> pl1 == pl2 False >>> pl1.equals(pl2) True For points, Python object equality works in some situations because of how simple points are spatially.  Once we move into polylines and polygons, and possibly multipoints, using Python object equality falls short. ... View more

How ArcPy is handling methods for testing spatial relations between geometric objects, e.g., equals, is consistent with the OGC Simple Feature Access - Part 1: Common Architecture use of Z and M coordinate values: 6.1.2.5 Use of Z and M coordinate values A Point value may include a z coordinate value.... A Point value may include an m coordinate value.... Observer methods returning Point values include z and m coordinate values when they are present. Spatial operations work in the "map geometry" of the data and will therefore not reflect z or m values in calculations (e.g., Equals, Length) or in generation of new geometry values (e.g., Buffer, ConvexHull, Intersection).  This is done by projecting the geometric objects onto the horizontal plane to obtain a "footprint" or "shadow" of the objects for the purposed of map calculations. In other words, it is possible to store and obtain z (and m) coordinate values but they are ignored in all other operations which are based on map geometries.  Implementations are free to include true 3D geometric operations, but should be consistent with ISO 19107. The operative statement for this discussion is really the last one, i.e., implementations are free to include true 3D [or 4D] geometric operations, but I will circle back to this in a minute. Other products that implement the OGC Simple Feature Access standards behave the same way as ArcPy in this instance.  Looking at Shapely, which is based on GEOS and JTS libraries, and grabbing a few example points from the original code snippets: >>> from shapely.geometry import Point >>> a = Point(1, 2, 3) >>> print a POINT Z (1 2 3) >>> c = Point(1, 2, 0) >>> print c POINT Z (1 2 0) >>> a == c False >>> a.equals(c) True Even looking at how MS SQL Server implements its Geometry data type: >>>Declare @a geometry, @c geometry >>>SET @a = 'POINT(1 2 3 4)' >>>PRINT @a.ToString() POINT (1 2 3 4) >>>SET @c = 'POINT(1 2 0 0)' >>>PRINT @c.ToString() POINT (1 2 0 0) >>>PRINT @a.STEquals(@c) 1 I think the behavior being seen in ArcPy is typical, some might even argue expected given the Simple Feature specifications/standards.  That said, Esri could implement OGC Feature Geometry/ ISO 19107 Geographic information -- Spatial schema in ArcPy, but even that would not necessarily address all of the issues here because Esri has created a 5D point model (x, Y, Z, M, ID). Has Esri failed or dropped the ball on this one?  It could be argued, but their existing implementation could also be defended on some levels.  The Simple Features specification sets the bar pretty low, but Esri also hasn't tried to raise it for themselves.  There is nothing in the specifications that prevents them from adding more robust functionality.  Think of how Microsoft handled the situation, i.e., they have OGC methods and Extended methods so they can remain compliant but add extra functionality. ... View more

In terms of null points, I completely agree with you that arcpy.Point() should not return the <Point (0.0, 0.0, #, #)> that it does.  I think Shapely/GEOS/JTS handles the situation much better: >>> from shapely.geometry import Point >>> null_pnt = Point() >>> print null_pnt GEOMETRYCOLLECTION EMPTY With an empty geometry collection, comparing a null/empty point to a valid regular point with 0, 0 coordinates behaves the way one would suspect: >>> real_pnt = Point(0, 0, 1) >>> print real_pnt POINT Z (0 0 1) >>> null_pnt.equals(real_pnt) False I think an empty geometry collection is the way to go; unfortunately, Esri hasn't implemented a geometry collection type with ArcPy.  Off the top, four different ideas come to mind how to handle the situation without a geometry collection: have arcpy.Point() throw an error saying the coordinates aren't valid implement an empty point concept and return an empty point return a None object return a valid point using the 0, 0 coordinates Out of the four options above, I think #4 (the existing condition) is the least desirable from a user/developer perspective. ... View more

I see you made progress with Xander Bakker​'s code, glad he was able to help and provide some additional learning links.  Beyond what has already been said, I will just point out that Python functions with 15+ parameters are quite rare, even 10 is unusual.  And, they are all required arguments!  In cases like this, usually the code can benefit from being restructured.  If all of these arguments have to be passed, maybe just building a list or dictionary and passing it as a single argument would make sense. ... View more

It is because you are returning with each if statement.  Once the first if condition is met, the code will build that label and exit the function. ... View more

In you want to start branching out into installing additional Python packages, I encourage you to read Installing Packages from the Python Packing User Guide.  ... View more

Try building the Extent object from scratch instead of reconstructing one you get from the DataFrame object.  If I understand what you are expecting for results, the following code works for me in the interactive Python window in ArcMap: custom_extent = {'XMin':-87.678956, 'XMax':-87.611652, 'YMin':41.867129, 'YMax':41.954932} newExtent = arcpy.Extent(**custom_extent) print data_frame.extent.XMin data_frame.spatialReference = arcpy.SpatialReference("WGS 1984") print data_frame.extent.XMin data_frame.extent = newExtent print data_frame.extent.XMin If you want to use a dictionary to store an extent's bounds, may I suggest using key names that correspond to function parameters.  That way, you can just pass the dictionary as keyword arguments to an Extent function. Update:  It dawned on me after posting the code above why your original code might not be working.  I believe it fails because you change the spatial reference of the data frame after getting the extent object from the data frame.  If you change the spatial reference first, it seems to work. custom_extent = {1:-87.678956, 2: -87.611652, 3:41.867129, 4:41.954932} print data_frame.extent.XMin data_frame.spatialReference = arcpy.SpatialReference("WGS 1984") print data_frame.extent.XMin newExtent = data_frame.extent newExtent.XMin, newExtent.XMax = custom_extent[1], custom_extent[2] newExtent.YMin, newExtent.YMax = custom_extent[3], custom_extent[4] data_frame.extent = newExtent print data_frame.extent.XMin ... View more

Since the OGC Simple Feature Access - Part 1: Common Architecture doesn't allow parts of a multipart polygon to share edges, what you are really after is either a geometry bag or geometry collection, neither of which are implemented in ArcPy.  Even if ArcPy implemented geometry bags, you still can't store them in a feature class because Esri doesn't have a geometry bag/collection type for feature classes. If you work outside of ArcGIS, in SQL Server for example, you can build and store a geometry collection.  Using your original example and a SQL Server workspace, you can see what I am speaking to: import arcpy # a square part = [[0,0],[0,1],[1,1],[1,0],[0,0]] mp = "" for x_shift in range(0,3):     for y_shift in range(0,5):         new_part = ", ".join(["{} {}".format(x+x_shift,y+y_shift) for [x,y] in part])         mp = "{}, POLYGON(({}))".format(mp, new_part) mp = mp[2:] sqlWS = #SQL Server or SQL Server Express SDE connection file out_name = "test_pol_mp" fc = arcpy.CreateFeatureclass_management(sqlWS, out_name, "POLYGON") shape = arcpy.Describe(fc).shapeFieldName sde_conn = arcpy.ArcSDESQLExecute(sqlWS) sql = ("INSERT INTO {} ({}) "       "VALUES (\'GEOMETRYCOLLECTION({})\')".format(out_name, shape, mp)) sde_conn.execute(sql) sql = ("SELECT [{}].STNumGeometries() AS partCount, "               "[{}].STNumPoints() AS pointCount "           "FROM {}".format(shape, shape, out_name)) partCount, pointCount = sde_conn.execute(sql)[0] print "Polygon partCount : {0}".format(partCount)  print "Polygon pointCount: {0}".format(pointCount) #... Polygon partCount : 15 #... Polygon pointCount: 75 Viewing the geometry collection in SQL Server Management Studio gives: Wheres trying to view it in ArcMap gives: ... View more

Darren Wiens you are correct in that it isn't an ArcPy limitation.  In fact, it isn't even an ArcGIS or Esri limitation, per se.  The quote from the Help link you provide is basically Esri paraphrasing an OGC standard without attribution or providing the 'why' to the user. The OGC Simple Feature Access - Part 1: Common Architecture document clearly states the "boundaries of any 2 Polygons that are elements of a MultiPolygon may not 'cross' and may touch at only a finite number of Points," i.e., parts of a multipart polygon can't share an edge.  MultiLineStrings can have shared lines or lines on top of each other.  The Esri Knowledge Base article FAQ:  Why are polygon features grouped into multi-line string features by the ArcSDE sdegroup command? speaks to this issue. In short, the standard doesn't allow for it.  Esri implements the standard, and their way of staying compliant is to dissolve parts of multipart polygons that share edges.  Microsoft implements the standard as well, but their way of handling the situation in SQL Server is to let the user create an invalid polygon and then let it err out later when someone tries to work with the geometry.  Same standard, two different products, and two different ways of coping with a user creating an invalid multipart polygon. ... View more

If the cursor isn't going to be used once and disposed of, using a Python with statement ensures the cursor is reset for its next use.  But then again, maybe having the cursor automatically reset isn't what someone wants.  Just depends on the situation and need. ... View more

True, for now and in the ArcPy realm, calling the search cursor's next method will get the same result.  I suggested using the built-in next method because of broader changes happening with Python outside of ArcPy.  For Python 3 after PEP 3114 was approved, it meant the next() iterator method was going away.  The Transition Plan for PEP 3114 covers two additional changes needed for moving to Python 3: Method definitions named next will be renamed to __next__ . Explicit calls to the next method will be replaced with calls to the built-in next function. For example, x.next() will become next(x) . The built-in next function was introduced in Python 2.6 to smooth the transition.  Since we know that Esri has made the leap to Python 3 with ArcGIS Pro, I suggested an approach using the built-in function.  For now, Esri's implementation of the ArcPy Data Access (arcpy.da) module in ArcGIS Pro still includes cursors having explicit next() methods, but I argue the current ArcGIS Pro implementation isn't very pythonic since the explicit next() methods don't add any special functionality beyond simply iterating. ... View more

Are you looking for only a single record/point per SRAddress?  I recall you mentioning there may be up to 10 types or collections.  Do you want a table that has [Type1, EWT1, ItemCount1, ...., Type10, EWT10, ItemCount10] columns?  A table with at least 30 columns, most of which will be empty? Or, are you wanting a separate table that has [SRAddress, Type, EWT, ItemCOunt] that you can use to relate back to a point feature class that has the point locations for each SRAddress? I guess I don't understand what you want the final data structure to look like in terms of how many tables and what those tables look like. ... View more

As Owen Earley‌ suggests, the Summary Statistics tool should work, and it is scriptable through Python. You could also use a data access (arcpy.da) search cursor approach, a bit more involved but likely more performant than Summary Statistics. # import functions from modules that are available but not commonly imported from collections import defaultdict from numpy import fromiter, dtype # sum PCE_VOLU_3 by NO stats = defaultdict(int) with arcpy.da.SearchCursor(input_fc, ['NO','PCE_VOLU_3']) as cur:     for k, v in cur:         stats += v # create iterable and populate numpy array  stats_iterable = ((k, v) for (k, v) in stats.iteritems())  tmp1 = fromiter(stats_iterable,                  dtype([('NO', 'i4'), ('SUM_PCE_', 'f8')])) # Dump numpy array to table  arcpy.da.NumPyArrayToTable(tmp1, out_table)  del tmp1  del stats ... View more

Once you sum the values, how do you want to store them?  Do you want a separate table?  Do you want a new column in this table where is shows the same sum for all of the records with the same NO? ... View more