康瑞部落

    网上有通过百度地图API获得特定城市的边界的程序，如：http://www.cnblogs.com/hzb2010/archive/2012/06/10/2544096.html能够批量下载县市级以上的城市边界；http://www.cnblogs.com/i-gps/archive/2012/05/18/2507941.html和http://www.cnblogs.com/i-gps/archive/2012/07/17/2595733.html以网页的形式下载区域的边界数据。以上两种办法使用的都是百度地图的数据源，区域的边界可能与真实的行政边界不一样。如：贵阳市南明区，其获得的边界如图1所示，但真实的南明区的范围如图2所示。

图1 已有方法获得的南明区边界

图2 贵阳市南明区的真实边界

    图2是从openstreetmap（https://www.openstreetmap.org/relation/2783470#map=11/26.5409/106.7765）中的查询结果。通过该网站，查询指定的区域，在Chrome开发者工具的network下，过滤“full”关键字，就能够获得边界数据了。下面是解析边界数据的代码：

static void ParseOpenStreetMap()
{
    string content = "";
    using (StreamReader sr = new StreamReader(@"下载下来的文件"))
    {
        content = sr.ReadToEnd();
    }
    XmlDocument doc = new XmlDocument();
    doc.LoadXml(content);
    Dictionary<string, List<double>> idToLatLng = new Dictionary<string, List<double>>();
    XmlNodeList dataNode = doc.SelectNodes("/osm/node");
    foreach (XmlNode node in dataNode)
    {
        string id = node.Attributes["id"].Value;
        double lat = double.Parse(node.Attributes["lat"].Value);
        double lng = double.Parse(node.Attributes["lon"].Value);
        idToLatLng[id] = new List<double>() { lat, lng };
    }

    dataNode = doc.SelectNodes("/osm/relation/member");
    List<string> wayOrder = new List<string>();
    foreach (XmlNode node in dataNode)
    {
        wayOrder.Add(node.Attributes["ref"].Value);
    }

    List<List<double>> latLngList = new List<List<double>>();
    HashSet<string> nodeIdSet = new HashSet<string>();
    int count = 0;
    foreach (var wayId in wayOrder)
    {
        dataNode = doc.SelectNodes("/osm/way[@id=" + wayId + "]/nd");

        List<List<double>> latLngListItem = new List<List<double>>();
        foreach (XmlNode node in dataNode)
        {
            string id = node.Attributes["ref"].Value;
            if (!nodeIdSet.Contains(id))
            {
                nodeIdSet.Add(id);
            }
            else
            {
                Console.WriteLine("dumplicated id in way/nd: " + id);
            }
            if (idToLatLng.ContainsKey(id))
            {
                latLngListItem.Add(idToLatLng[id]);
            }
            else
            {
                Console.WriteLine("no id: " + id);
            }
        }
        if (count == 1)    //第一个放下来的可能是反的
        {
            if (GetDistance(latLngList.Last()[0], latLngList.Last()[1], latLngListItem.First()[0], latLngListItem.First()[1]) > 1
                &&
                GetDistance(latLngList.Last()[0], latLngList.Last()[1], latLngListItem.Last()[0], latLngListItem.Last()[1]) > 1)
            {
                latLngList.Reverse();
            }
        }
        if (latLngList.Count == 0 || GetDistance(latLngList.Last()[0], latLngList.Last()[1], latLngListItem.First()[0], latLngListItem.First()[1]) <= 1)
        {
            latLngList.AddRange(latLngListItem);
        }
        else
        {
            latLngListItem.Reverse();
            latLngList.AddRange(latLngListItem);
        }
        count++;
    }
                        
    using (StreamWriter sw = new StreamWriter(@"解析后的输出文件"))
    {
        foreach (var l in latLngList)
        {
            sw.Write("[" + l[0] + "," + l[1] + "],");
        }
    }
    Console.WriteLine(latLngList.Count);
}
private const double EARTH_RADIUS = 6378.137; //地球半径
private static double rad(double d)
{
    return d * Math.PI / 180.0;
}

public static double GetDistance(double lat1, double lng1, double lat2, double lng2)
{
    double radLat1 = rad(lat1);
    double radLat2 = rad(lat2);
    double a = radLat1 - radLat2;
    double b = rad(lng1) - rad(lng2);
    double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) +
    Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Pow(Math.Sin(b / 2), 2)));
    s = s * EARTH_RADIUS;
    s = Math.Round(s * 10000) / 10000;
    return s;    //单位：千米
}