6

Basically I have set of points [lat-long pairs] which define different polygons on a map. Now I have a coordinate [lat-long] from which I want to find the nearest polygon from that point.

  • Do you specifically want this with google-maps?, you could this right away with ArcMap or QGIS – Daniel Jul 10 at 9:38
7

There are two approaches. One is brute-force: write code to find the distance between a point and a polygon, and then iterate over all polygons. The second involves building efficient data structures, an arduous process.

For the first approach, you might look at this Wolfram Demonstration Project for computing the distance from a point to a polygon. Link here; image below. Note that sometimes the min distance is achieved at a point interior to an edge of the polygon.
           Min distance to polygon

  • +1 for answering question... Can these be achieved using PHP? – Rajeev Vyas Jan 2 '13 at 12:41
  • Of course, Yes, assuming you mean the php scripting language. It is singularly inappropriate for this type of computation, but one could do it. I would instead use a language more oriented toward numerical computation. – Joseph O'Rourke Jan 3 '13 at 0:27
  • @Joseph O'Rourke: Care to elaborate about "building efficient data structures" in this context ? – Stefan Steiger Jul 10 at 9:44
  • @StefanSteiger: I had in mind computing the Voronoi diagram of edges and vertices, after which queries are efficient (sublinear). – Joseph O'Rourke Jul 11 at 10:39
0

Necromancing.

So, assuming you have a point (geocoded coordinates of building address).
Now you want to get the building polygon.

There are several problems with this. First, nominatim/OSM is crap in geocoding, so you have to use the google geocoding API, which is not free anymore. Alternatively you could use the yandex geocoding API.

Second, google has no API that gives you the building polygons. So you'll either need to vectorize google-maps-tiles, or you need to take the polygons from the XML-API of OpenStreetMap.

For vectorizing google-maps polygons, I'll relegate you to this method using photoshop, and this method using python.

Since that's obviously crap, I'm going to assume you're using the OSM polygon API.

Which leads us to to problem number 3:
Google-Geocoding API may only be used in conjuction with Google-Maps.
So if you're using it with OSM, that's against the terms of service, strictly speaking (aka illegal).
Same shenanigan is true for the Yandex-API.
Besides, the yandex-API gives you points on the street, not a point inside the building-polygon, and not the barycenter of the building-polygon either (last one also true for google-maps).

Also, note that OSM-polygons do not necessarely match with google-polygons (displacement / off by x-meters errors).

So, assuming you have been able to get a hold of a manually/legally obtained WGS-84 point (building address coordinate), that you may use with OSM, you can do this:

Get all OSM-data inside a certain bounding-box:

const string OSM_API_VERSION = "0.6";
// string url = "https://www.openstreetmap.org/api/0.6/map?bbox=8.626273870468141,47.69679769756054,8.636573553085329,47.700530864557194&no_cache=1562588642802";
string url = "https://www.openstreetmap.org/api/" + OSM_API_VERSION + "/map?bbox=" + bounds.ToBBoxString();

string xml = null;
using (System.Net.WebClient wc = new System.Net.WebClient())
{
    xml = wc.DownloadString(url);
} // End Using wc

where

public string ToBBoxString()
{
    decimal[] a = new decimal[]{this.West, this.South, this.East, this.North};

    return a.Join(",");
}

Note: The bounding-box area my not exceed 0.25 old-degrees squared !

decimal area = bounds.BoundsArea;
if (area > 0.25m)
{
    System.Console.WriteLine("The maximum bbox size is 0.25, and your request was too large.\nEither request a smaller area, or use planet.osm.");
    return;
} // End if (area > 0.25m) 

where

public decimal BoundsArea
{
    get
    {
        LatLng nw = this.NorthWest;
        LatLng se = this.SouthEast;
        // https://github.com/openstreetmap/cgimap/blob/master/src/bbox.cpp

        decimal maxLng = System.Math.Max(nw.lng, se.lng);
        decimal maxLat = System.Math.Max(nw.lat, se.lat);

        decimal minLng = System.Math.Min(nw.lng, se.lng);
        decimal minLat = System.Math.Min(nw.lat, se.lat);
        decimal area = (maxLng - minLng) * (maxLat - minLat);

        return area;
    }

} // End Property BoundsArea 

To construct the bounding box, take your WGS-84 point, and make it a 1000 meters bounding box (that will get you everything within a radius of 500 meters).

Like this:

public static LatLngBounds FromPoint(LatLng point, decimal sizeInMeters)
{
    decimal latAccuracy = 180.0m * sizeInMeters / 40075017m;
    decimal lngAccuracy = latAccuracy / (decimal)System.Math.Cos((System.Math.PI / 180.0d) * (double)point.lat);

    //           N
    //          180
    // (W) -180     +180 (E)
    //         -180
    //           S

    // https://github.com/Leaflet/Leaflet/blob/master/src/geo/LatLng.js
    // https://github.com/Leaflet/Leaflet/blob/master/src/geo/LatLngBounds.js
    // constructor(southWest: LatLngExpression, northEast: LatLngExpression);
    // a = [point.lat - latAccuracy, point.lng - lngAccuracy],
    // b = [point.lat + latAccuracy, point.lng + lngAccuracy]
    // new LatLngBounds(a, b); 

    decimal south = point.lat - latAccuracy;
    decimal west = point.lng - lngAccuracy;
    decimal north = point.lat + latAccuracy;
    decimal east = point.lng + lngAccuracy;


    // https://en.wikipedia.org/wiki/Ellipse
    // https://en.wikipedia.org/wiki/Latitude
    // https://en.wikipedia.org/wiki/Longitude
    // https://de.wikipedia.org/wiki/Wendekreis_(Breitenkreis)

    return new LatLngBounds(north, south, east, west);
}

Now that you have all the nearby polygon's data, you'll need to construct the polygons from it, and then measure for every polygon side and every polygon edge-point, the distance to your WGS-84 point. The polygon you want is the nearest polygon to your point (assuming it exists - otherwise it takes the nearby-est), so it will be the polygon with the smallest distance to your WGS-84 point.

This goes like this:

System.Xml.XmlDocument doc = new System.Xml.XmlDocument();
doc.LoadXml(xml);

System.Xml.XmlNodeList nodes = doc.SelectNodes("//node");


System.Collections.Generic.Dictionary<string, LatLng> nodeDictionary =
    new System.Collections.Generic.Dictionary<string, LatLng>(System.StringComparer.InvariantCultureIgnoreCase);

System.Collections.Generic.Dictionary<string, LatLng[]> buildingDictionary =
    new System.Collections.Generic.Dictionary<string, LatLng[]>(System.StringComparer.InvariantCultureIgnoreCase);


foreach (System.Xml.XmlElement node in nodes)
{
    string id = node.GetAttribute("id");
    string nodeLat = node.GetAttribute("lat");
    string nodeLong = node.GetAttribute("lon");

    decimal dlat = 0;
    decimal dlong = 0;
    decimal.TryParse(nodeLat, out dlat);
    decimal.TryParse(nodeLong, out dlong);

    nodeDictionary[id] = new LatLng(dlat, dlong);
} // Next node 

// https://stackoverflow.com/questions/1457638/xpath-get-nodes-where-child-node-contains-an-attribute
// querySelectorAll('way tag[k="building"]')
System.Xml.XmlNodeList buildings = doc.SelectNodes("//way[tag/@k=\"building\"]");
foreach (System.Xml.XmlElement building in buildings)
{
    System.Collections.Generic.List<LatLng> lsPolygonPoints = new System.Collections.Generic.List<LatLng>();

    System.Xml.XmlNodeList buildingNodes = building.SelectNodes("./nd");
    foreach (System.Xml.XmlElement buildingNode in buildingNodes)
    {
        string reff = buildingNode.GetAttribute("ref");
        lsPolygonPoints.Add(nodeDictionary[reff]);
    } // Next buildingNode 


    LatLng[] polygon = toCounterClockWise(lsPolygonPoints.ToArray());
    string id = building.GetAttribute("id");
    buildingDictionary[id] = polygon;

    string sqlPolygon = CreateSqlPolygon(polygon);
    System.Console.WriteLine(sqlPolygon);
} // Next building 

System.Console.WriteLine(buildingDictionary);


System.Collections.Generic.Dictionary<string, Polygon> buildingPolygonDictionary =
    new System.Collections.Generic.Dictionary<string, Polygon>(System.StringComparer.InvariantCultureIgnoreCase);

foreach (System.Collections.Generic.KeyValuePair<string, LatLng[]> kvp in buildingDictionary)
{
    buildingPolygonDictionary[kvp.Key] = new Polygon(kvp.Value);
} // Next kvp 

System.Console.WriteLine(buildingPolygonDictionary);



decimal? min = null;
string uid = null;

foreach (System.Collections.Generic.KeyValuePair<string, Polygon> kvp in buildingPolygonDictionary)
{
    decimal minDist = kvp.Value.GetMinimumDistance(geoPoint);

    if (min.HasValue)
    {
        if (minDist < min.Value)
        {
            min = minDist;
            uid = kvp.Key;
        }
    }
    else
    {
        uid = kvp.Key;
        min = minDist;
    }

} // Next kvp 

Polygon p = buildingPolygonDictionary[uid];
System.Console.WriteLine(uid); // 218003784
System.Console.WriteLine(p);


LatLng[] pts = p.ToLatLngPoints();
string sql = CreateSqlPolygon(pts);
System.Console.WriteLine(sql);

And now you have the nearest polygon.
Of course, all the really simple vector geometry is now hidden inside

DecimalVector2.GetPointVerticalIntersection 

which is called from

kvp.Value.GetMinimumDistance(geoPoint);

which is not helping you. I could now leave you like that, but since I'm not that cruel, I'll instead supply you with the class DecimalVector2 and Polygon, for reference, if you need to implement it in another programming language.

Note: If you need to save the polygon in the database, know that MS-SQL requires the polygon's points to be of CLOCKWISE sequence, while to render them on leaflet/OSM, they need to be in COUNTER-clockwise sequence. To that end, the class Polygon has methods ToClockWiseLatLngPoints and ToCounterClockWiseLatLngPoints.

Also, the polygon's points need to be closed, aka

point[0] = point[N]

or

public bool IsClosed
{
    get
    {
        decimal x0 = this.Points[0].X;
        decimal y0 = this.Points[0].Y;

        decimal xN = this.Points[this.Points.Count - 1].X;
        decimal yN = this.Points[this.Points.Count - 1].Y;

        if (x0 == xN && y0 == yN)
            return true;

        return false;
    }
}



public void Close()
{
    Points.Add(new DecimalVector2(this.Points[0].X, this.Points[0].Y));
}

public void UnClose()
{
    Points.RemoveAt(Points.Count - 1);
}


public void EnsureClosed()
{
    if (!this.IsClosed)
        this.Close();
}

Also, not that you're just comparing distances here.
That means you only need relative comparisons.
So in every instance where sqrt(x^2 +y^2) is used (aka vector-length), you could omit the SQRT to make the program faster/more-efficient.

Currently, with square-root and 1000 meters bounding-box, it takes about 1.5 seconds per building (in debug mode).

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.