【原】数据清洗告一段落了

老马的程序人生 2020-08-17

展开全文

今天8:00开始流水判线性代数试卷，直到17:30，经过2天半的时间，终于批改完毕，剩下的就是核对总分提交成绩了。所以数据清洗的任务只能放到晚上来进行，这次主要是把数据从大地坐标系（经度、纬度）转化到平面坐标系（X，Y）。

经度是指某点与两极的连线与0度经线所在平面的夹角，国际上规定以通过英国伦敦近郊的格林尼治天文台旧址的经线作为计算经度的起点，即经度零度零分零秒，也称“本初子午线”。它东面的为东经，记为E，共180度，西面的为西经，记为W，共180度。
纬度是指某点与地球球心的连线和地球赤道面所成的线面角，其数值在0至90度之间。位于赤道以北的纬度称为北纬，记为N；位于赤道以南的纬度称为南纬，记为S。

通常情况下我们利用 Baidu API 函数得到的坐标系是 DB09 坐标系，在进行投影变换之前，要转换为国际通用的 WGS84坐标系。

WGS84：为一种大地坐标系，也是目前广泛使用的GPS全球卫星定位系统使用的坐标系。
GCJ02：又称火星坐标系，是由中国国家测绘局制定的地理坐标系统，是由WGS84加密后得到的坐标系。
BD09：为百度坐标系，在GCJ02坐标系基础上再次加密。其中bd09ll表示百度经纬度坐标，bd09mc表示百度墨卡托米制坐标。

坐标系转换的代码如下：

public class Coordtransform
{
    private const double XPi = 3.14159265358979324 * 3000.0 / 180.0;
    private const double Pi = 3.1415926535897932384626;// π
    private const double A = 6378245.0;// 长半轴
    private const double Ee = 0.00669342162296594323;// 偏心率平方

    // 百度坐标系 (BD-09) 与 火星坐标系 (GCJ-02)的转换
    public static double[] Bd09ToGcj02(double bdLon, double bdLat)
    {
        double x = bdLon - 0.0065;
        double y = bdLat - 0.006;
        double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * XPi);
        double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * XPi);
        double ggLng = z * Math.Cos(theta);
        double ggLat = z * Math.Sin(theta);

        return new[] { ggLng, ggLat };
    }

    // 判断是否在国内，不在国内则不做偏移纬度3.86~53.55,经度73.66~135.05
    public static bool Out_of_china(double lng, double lat)
    {
        return !(lng > 73.66 && lng < 135.05 && lat > 3.86 && lat < 53.55); ;
    }

    public static double Transformlat(double lng, double lat)
    {
        double ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * Math.Sqrt(Math.Abs(lng));
        ret += (20.0 * Math.Sin(6.0 * lng * Pi) + 20.0 * Math.Sin(2.0 * lng * Pi)) * 2.0 / 3.0;
        ret += (20.0 * Math.Sin(lat * Pi) + 40.0 * Math.Sin(lat / 3.0 * Pi)) * 2.0 / 3.0;
        ret += (160.0 * Math.Sin(lat / 12.0 * Pi) + 320 * Math.Sin(lat * Pi / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    public static double Transformlng(double lng, double lat)
    {
        double ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.Sqrt(Math.Abs(lng));
        ret += (20.0 * Math.Sin(6.0 * lng * Pi) + 20.0 * Math.Sin(2.0 * lng * Pi)) * 2.0 / 3.0;
        ret += (20.0 * Math.Sin(lng * Pi) + 40.0 * Math.Sin(lng / 3.0 * Pi)) * 2.0 / 3.0;
        ret += (150.0 * Math.Sin(lng / 12.0 * Pi) + 300.0 * Math.Sin(lng / 30.0 * Pi)) * 2.0 / 3.0;
        return ret;
    }


    // GCJ02 转换为 WGS84
    public static double[] Gcj02ToWgs84(double lng, double lat)
    {
        if (Out_of_china(lng, lat))
        {
            return new[] { lng, lat };
        }
        else
        {
            double dlat = Transformlat(lng - 105.0, lat - 35.0);
            double dlng = Transformlng(lng - 105.0, lat - 35.0);
            double radlat = lat / 180.0 * Pi;
            double magic = Math.Sin(radlat);
            magic = 1 - Ee * magic * magic;
            double sqrtmagic = Math.Sqrt(magic);
            dlat = (dlat * 180.0) / ((A * (1 - Ee)) / (magic * sqrtmagic) * Pi);
            dlng = (dlng * 180.0) / (A / sqrtmagic * Math.Cos(radlat) * Pi);
            double mglat = lat + dlat;
            double mglng = lng + dlng;

            return new[] { lng * 2 - mglng, lat * 2 - mglat };
        }
    }

    // 百度坐标系->WGS84坐标系
    public static double[] Bd09ToWgs84(double bdLon, double bdLat)
    {
        double[] gcj02 = Bd09ToGcj02(bdLon, bdLat);

        return Gcj02ToWgs84(gcj02[0], gcj02[1]);
    }
}