9

This is ArcGIS specific.

I have 2 point shapefiles A and B, the first one (A) is a single point that contains a lat long, the second (B) is a multitude of points (over 12k) that each contain their lat and long. What I am trying to do is to automate the selection of 75% of shapefile B's points based on the distance from shapefile A. In other words I want to select the nearest 75% of shapefile B points to shapefile A's one point.

5
  • Is a programmatic solution acceptable? May 5, 2011 at 23:16
  • BTW, I requested that Esri allow the Shapefield to be used in a custom ITableSortCallback but was told there was no justification for this. This use case shows otherwise. May 5, 2011 at 23:19
  • @Kirk Kuykendall Yes a programmatic solution would actually be preferred as this is a process I will have to repeat over 1k times. I have roughly 1200 separate points and each of those points has another shapefile with on average 12k points around it. I need to figure out a way to easily select the nearest 75% of surrounding points for them all. Doing it manually is just out of the question.
    – Furlong
    May 6, 2011 at 15:48
  • Perhaps this comment is outside the proper scope of a comment, but when and why would such an analysis be useful? This is for my own elucidation; forgive my slowness.
    – Nathanus
    May 6, 2011 at 18:41
  • 1
    Consider using statistical software. If you were to merge all 1200 shapefiles, creating a source id field during the merge, you could join the corresponding central point coordinates to that and compute all 1200 * 12k = 14.4M distances. What you need then is a list of the 75th percentiles of distance by source id: that would take about ten seconds with Stata (commercial) or R (open source). (If you use ArcGIS for this, let us know how much time it takes for the computation. :-)
    – whuber
    May 6, 2011 at 20:20

4 Answers 4

5

You could make a Multiple Ring Buffer on shapefile A, and then do a spatial join of the Buffer to shapefile B. When you do a spatial join of polygons and points, you get a count of the number of points in each polygon in the attribute table of the join. Then, by examining the total number of points within the buffers, you can get within 75% of the points in shapefile B.

A slightly different approach would be to script this in Python, and check for the 75% in a loop, but if it is a once off calculation you may not need that.

8
  • 4
    It would be simpler, faster, and more accurate to perform a spatial join of A to B, compute the third quartile of the resulting [distance] field, and select all records less than that distance.
    – whuber
    May 6, 2011 at 4:32
  • I did not know it was possible to spatially join points! I agree, this would be a much better way of doing it.
    – djq
    May 6, 2011 at 10:54
  • @Andy On the contrary, the join is a nearest-point relation. It is not based on any tabulated attributes at all. Also, in Arc* software (going back to ArcView 2), the distance is automatically computed as a result of the join.
    – whuber
    May 6, 2011 at 12:33
  • 1
    @whuber, I know! Hence the retracted (deleted statement!) I assume you could do it by attribute table joins (and calculating the distance yourself) but that would be unnecessary given the context. I guess the point I would like to reiterate is that it is simply calculating the distance between 1 point, no looping or buffers or iterative procedures are necessary.
    – Andy W
    May 6, 2011 at 13:36
  • 1
    @Furlong If you read the example of Spatial Join: help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… you can get an idea how to run this in python. Then, it's a matter of running through the attribute table and choosing values that match your criteria
    – djq
    May 6, 2011 at 16:12
4

For 1200 points (or even up to say 12M points?) I'd just put them into memory as a Generic Collection - in this case a SortedList of lists. This could be simplified by just skipping points when you run into a situation with multiple points that are same distance from the origin point. Also, for performance, consider using a hashtable instead of a SortedList, and sorting once after inserting all distances. That would take a few more lines of code though(?).

I didn't have time to test this, but this c# might get you started:

private void SelectNTile(string layer1, string layer2, double nTile)
{
    var fLayer1 = FindLayer(ArcMap.Document.FocusMap, "LayerWithLotsofPoints");
    var fLayer2 = FindLayer(ArcMap.Document.FocusMap, "LayerWithOneSelectedPoint");
    IFeature feat = GetSingleFeature(fLayer2);
    var distList = MakeDistList(fLayer1.FeatureClass,(IPoint)feat.ShapeCopy);
    // assume not many points exactly same distance
    var nRecs = (int)(distList.Count * nTile); // nTile would be 0.75 for 75%
    var Oids = new List<int>();
    foreach (KeyValuePair<double, List<int>> kvp in distList)
    {
        Oids.AddRange(kvp.Value);
        if (Oids.Count > nRecs)
            break;
    }
    var fSel = fLayer1 as IFeatureSelection;
    var OidArray = Oids.ToArray();
    fSel.SelectionSet.AddList(Oids.Count, ref OidArray[0]);                
}

private SortedList<double, List<int>> MakeDistList(IFeatureClass fc, IPoint pnt)
{
    var outList = new SortedList<double, List<int>>();
    var proxOp = pnt as IProximityOperator;
    IFeatureCursor fCur = null;
    try
    {
        fCur = fc.Search(null, true); // recycling is faster, we just need OIDs
        IFeature feat;
        while ((feat = fCur.NextFeature()) != null)
        {
            double dist = proxOp.ReturnDistance(feat.Shape);
            if (!outList.ContainsKey(dist))
                outList.Add(dist, new List<int> { feat.OID });
            else
                outList[dist].Add(feat.OID);  // this should rarely happen
        }
    }
    catch
    {
        throw;
    }
    finally
    {
        if (fCur != null)
            System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
    }
    return outList;
}
private IFeature GetSingleFeature(IFeatureLayer fLayer)
{
    var fSel = fLayer as IFeatureSelection;
    if (fSel.SelectionSet.Count != 1)
        throw new Exception("select one feature in " + fLayer.Name + " first");
    var enumIDs = fSel.SelectionSet.IDs;
    enumIDs.Reset();
    IFeature feat = fLayer.FeatureClass.GetFeature(enumIDs.Next());
    return feat;
}
private IFeatureLayer FindLayer(IMap map, string name)
{
    throw new NotImplementedException();
}
4

A Python geoprocessing script is an obvious choice:

  1. Use the Point Distance tool to calculate the distance from your features in feature class B (the "Input Features" parameter of the tool) to the point in feature class A (the "Near Features" parameter of the tool).
  2. Sort the table by the calculated distance.
  3. Select the first 75% of the objectids in the output table (the "Input_FID" column) and use those to make your selection from the original features in feature class B.
2

I had this issue a few years back. I found it easier to keep the data as 'flat data', looping through all of the data and manually computing the distance, then taking the top 75% (I actually kept the top 10%). I then did the same in ArcIMS using their distance calculations and it took a lot longer.

Buffering is a huge overhead, yet maths calculations are a 'puters forte. If you buffer 12k points, I think you're going to have performace issues.

1
  • I [@Mapperz] deleted the comments - mod tools guidelines flagged this post because it degraded into pointless bickering...
    – Mapperz
    May 20, 2011 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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