1

I want to find the geodesic distance between the closest hydrography feature to each of my lat/long coordinates.

I have hydrography data from New York state in a shapefile. I have many points, 10^10 number of coordinates in WGS1984 datum. There's about 20,000 features in my hydrography data which was originally in NAD1983.

I'm trying to find a way to expedite my code. This is what I've tried:

Near function in ArcMap: This worked reasonably well when I used roads (about 10,000 features). It still takes about a week but is taking appears to be taking forever with hydrography.

Potential Thoughts:

I've converted all of my hydrography features to WGS 1984 to work with it. Should I consider converting both files into planar coordinates? I'm not sure I understand well enough whether converting from unprojected to projected would cause any potential problems.

I've tried methods in both ArcMap or R, I'm open to either and just want the fastest way to process all of these points. I've written a second question in the event somebody has a better method in R.

Expedite Near function in R for 10^9 coordinates?

  • Ten billion points? 10k vs 20k lines, with how many vertices? Why shapefile, and not a modern format like file geodatabase? How large an area? What search distance? Did you dice the features to reduce false positives? Have you tried in_memory storage? – Vince Dec 11 '19 at 15:09
  • I would recommend re-projecting to a projected coordinate system, import the layers into a file geodatabase, and add spatial index to each table before running near in ArcGIS. – artwork21 Dec 11 '19 at 15:28
  • 2
    By making this question both ArcGIS and R, you make it two questions, in violation of the One question per Question policy. Since the comments are focused on the ArcGIS solution, I suggestion you focus the question on just ArcGIS, and create a second question for how to optimize R. – Vince Dec 11 '19 at 15:45
  • Hi Vince, Thanks for your comment--I didn't think that would be considered one question. I'll go ahead and split it into two. Apologies, I should have put 10^9 (about 200 million points). I'm not quite sure that I'm following correctly--the hydrography is in a shapefile and the lat/long are split into text files, which I've added + saved as shapefiles, and then run Near function with each one. If I understand correctly, are you suggesting converting all my points into a geodatabase instead? I'm searching across all of New York State. Could you explain "dice the features" and "in_memory"? – Tammy Dec 11 '19 at 16:31
  • 1
    The computational complexity of a Near operation is O(number_of_vertices). A roads feature class with 10,000 features could have 20,000 vertices or 200 billion. When comparing the expected difference between a 10k feature shapefile and a 20k feature shapefile, it is necessary to report the difference in the number of vertices, not the number of features. While not necessary to approach a solution, reporting the difference in vertex counts might account for the difference in execution time. To calculate, you'd need to iterate a DA SearchCursor, summing the shape pointcount property. – Vince Dec 11 '19 at 21:00
3

I generated a bunch of data, did some queries and found some interesting results.

First, I started with 1137 hydrographic features (27,232 vertices) for New York (geographic CS, NAD83 datum), in shapefile hydro. Then I densified the vertices at a 100 meter interval to kick the vertex count up to 187,525 in shapefile hydro_100m.

Next I generated hexagon tessellations at 0.001 and 0.0001 degree separation, selected for features within NY state (1:25m scale), and generated centroid points for those hexagons, yielding point shapefiles points_e03 and points_e04, with 14,414 and 140,409 features, respectively.

Then I copied the hydro, hydro_100m, points_e03, and points_e04 shapefiles to file geodatabase, and hydro, hydro_100m to the in_memory workspace.

The runtimes in seconds for base shapefiles against shapefiles using GEODESIC distance calculation was:

                         points_e03.shp  points_e04.shp 
    hydro.shp               9.46            90
    hydro_100m.shp         17.09           163

Then I used file geodatabase hydrography, but that didn't make an improvement:

                         points_e03.shp  points_e04.shp 
    hydro                   9.31            90
    hydro_100m             17.61           167

More surprisingly, even in_memory hydrography didn't improve performance:

                         points_e03.shp  points_e04.shp 
    in_memory/hydro         9.64            92
    in_memory/hydro_100m   17.55           167

When I changed the points to file geodatabase, there was a slight performance improvement (in the smaller point table, at least):

                         points_e03      points_e04 
    hydro.shp               8.87            92
    hydro_100m.shp         17.02           165

                         points_e03      points_e04 
    hydro                   8.54            90
    hydro_100m             16.47           165

                         points_e03      points_e04 
    in_memory/hydro         9.12            91
    in_memory/hydro_100m   17.10           165

Then I tried customizing the FGDB spatial reference to limit the precision to 1.0e-07 degrees:

                         opt_points_e03  opt_points_e04 
    opt_hydro               8.74            84
    opt_hydro_100m         16.63           158

                         opt_points_e03  opt_points_e04 
    in_memory/ohydro        8.58            82
    in_memory/ohydro_100m  16.22           155

So then I projected into the USGS Albers Equal Area the contiguous US (with 1 centimeter precision), and tried again using Cartesian distance calculations (PLANAR option):

                         apoints_e03     apoints_e04    
    ahydro.shp              1.14             9.78
    ahydro_100m.shp         1.56            11.38

                         apoints_e03     apoints_e04    
    ahydro                  1.11             9.56
    ahydro_100m             1.55            11.11

                         apoints_e03     apoints_e04    
    in_memory/ahydro        1.11             9.73
    in_memory/ahydro_100m   1.53            11.22

Yowza! Now we're talking!

Just for grins, I tried using Dice on the Albers hydro feature classes, to gauge impact of smaller features (with smaller bounding rectangle, which feeds into index performance). The ahydro layer only started with 23 mean vertices, so I diced it to 12, and the ahydro_100m started with ~165 mean vertices, and I diced it to 30 (the actual new means were 9.85 and 27.62, respectively, which are pretty small). Then another pass:

                         apoints_e03     apoints_e04    
    dhydro.shp              1.19             9.74
    dhydro_100m.shp         1.72            11.36

                         apoints_e03     apoints_e04    
    dhydro                  1.15             9.91
    dhydro_100m             1.70            12.21

                         apoints_e03     apoints_e04    
    in_memory/dhydro        1.17            10.07
    in_memory/dhydro_100m   1.56            11.47

(Okay, so you can't win 'em all, but for contours and other features that do spread across large envelopes, the Dice tool can be a benefit.)

And, since my points were created by a regular generator, let's try one more time, with the same points in random physical order (sorted by a random field, which was then dropped). This should show if spatial fragmentation of the point source has an impact on Near calculation:

                         rpoints_e03     rpoints_e04    
    ahydro.shp              1.19            10.35
    ahydro_100m.shp         1.62            13.19

                         rpoints_e03     rpoints_e04    
    ahydro                  1.24            11.09
    ahydro_100m             1.63            12.80

                         rpoints_e03     rpoints_e04    
    in_memory/ahydro        1.19            10.53
    in_memory/ahydro_100m   1.66            12.77

So, lessons learned:

  • It appears that the Near command may be placing the linear features in a cache, so that the Near performance is optimized, no matter the source (see Caveat).
  • Reprojecting had a tremendous impact (order of magnitude!)
  • Shapefile doesn't hurt performance as much as it might, though if there were lots of string attributes, then the raw I/O of rewriting the result rows might have had an impact
  • Using optimized file geodatabase spatial references had a slight (measurable) impact, but was still far less than elimination of GEODESIC
  • The Dice utility didn't help this particular dataset
  • Spatial fragmentation rears its ugly head on all large table operations, even as small as 140k rows, so processing the points to be spatially organized is likely to help in massive table operations

Caveat: My work laptop is not a slacker -- It has 4 core x 2.9Ghz Intel i7 CPU, 16GB RAM, and 2x1000GB SSD, so it's possible the sub-millisecond solid-state disk seek pushed the differences between disk and in_memory out of the measurable range. If you're not using an SSD for processing a 200 million row table, you might have to add a zero or two before the decimal place of the expected runtime.

FWIW: The code I used to calculate mean vertex counts was just pasted into the Python window of ArcMap:

feats = 0
verts = 0
with arcpy.da.SearchCursor("dhydro_100m",['shape@']) as cursor:
    for row in cursor:
        feats += 1
        verts += row[0].pointCount
print("feats = {:d}, verts = {:d}, mean = {:.2f}".format(
    feats,verts, float(verts)/float(feats)))

  • Wow, can't even thank you for your time & didactic answer---this was a large mountain to start climbing and makes this more manageable. To check on your previous comment and my understanding, using projected coordinates will essentially ignore the distortion from the curvature of the earth. Albers is best around middle latitudes (reading about it, it sounds like I shouldn't use it past Lat 35 deg). I would have thought NYS is a large enough area to use geodesic distances. Could you shed some insight on this? I just want to understand the caveats for what sounds like the only feasible option. – Tammy Dec 12 '19 at 22:19
  • I just used a nominal projection to measure performance difference. You should use an appropriate custom projection, one which manages the distortion. – Vince Dec 13 '19 at 0:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.