技術文章 » 經緯度座標為中心點生成米距離長度半徑的圓形面，含java js原始碼+線上繪製，程式碼簡單零依賴

經緯度座標為中心點生成米距離長度半徑的圓形面，含java js原始碼+線上繪製，程式碼簡單零依賴

2023-06-08 12:01:18

目錄

前些時間在更新我的座標邊界查詢工具的時候，需要用到經緯度座標點的距離計算，和以座標點為中心生成一個指定距離為半徑的圓，搜了一下沒有找到現成簡單又合適的程式碼，於是把自己壓箱底的程式碼翻出來了，簡化完善了一下，嘿，程式碼量也不大，還挺好用。

本方法是通過計算得到圓上的多個座標點，來得到的一個近似的圓形面，只要座標點夠多，這個圓就能足夠圓；有了這些座標點就很容易表示成不同的格式，比如：GeoJSON文字、WKT文字、Geometry範例。

源自 座標邊界查詢工具 開源庫：https://github.com/xiangyuecn/AreaCity-Query-Geometry （github可以換成gitee），高效能的座標資料、邊界資料查詢工具，Java開源程式、帶http查詢介面，記憶體佔用低（1秒可查1萬個以上座標對應的城市資訊）。

java版原始碼

public static void main(String[] args) {
	//計算天壇到天安門的距離
	System.out.println(Distance(116.410622, 39.881773, 116.397476, 39.908647));
	//生成天壇1公里範圍的圓形面
	System.out.println(CreateSimpleCircleWKT(116.410622, 39.881773, 1000, 24));
}

/** 計算兩個座標的距離，單位米 **/
static public double Distance(double lng1, double lat1, double lng2, double lat2) {
	//採用Haversine formula演演算法，高德地圖的js計算程式碼，比較簡潔 https://www.cnblogs.com/ggz19/p/7551088.html
	double d=Math.PI/180;
	double f=lat1*d, h=lat2*d;
	double i=lng2*d - lng1*d;
	double e=(1 - Math.cos(h - f) + (1 - Math.cos(i)) * Math.cos(f) * Math.cos(h)) / 2;
	return 2 * 6378137 * Math.asin(Math.sqrt(e));
}

/** 以座標點為中心，簡單粗略的建立一個指定半徑的圓，半徑單位米，pointCount為構建圓的座標點數（比如24個點，點越多越圓，最少3個點），返回構成圓的座標點陣列 **/
static public double[][] CreateSimpleCircle(double lng, double lat, double radius, int pointCount) {
	//球面座標不會算，轉換成三角座標簡單點，經度代表值大約：0.01≈1km 0.1≈10km 1≈100km 10≈1000km
	double km=radius/1000;
	double a=km<5?0.01 :km<50?0.1 :km<500?1 :10;
	double b=Distance(lng, lat, lng+a, lat);
	double c=Distance(lng, lat, lng, lat+a);
	double rb=radius/b*a;
	double rc=radius/c*a;
	double[][] arr=new double[pointCount+1][];
	double n=0,step=360.0/pointCount,N=360-step/2; //注意浮點數±0.000000001的差異
	for(int i=0;n<N;i++,n+=step){
		double x=lng+rb*Math.cos(n*Math.PI/180);
		double y=lat+rc*Math.sin(n*Math.PI/180);
		arr[i]=new double[] { x, y };
	}
	arr[pointCount]=new double[] { arr[0][0], arr[0][1] }; //閉環
	return arr;
}

/**
以座標點為中心，簡單粗略的建立一個指定半徑的圓，半徑單位米，pointCount為構建圓的座標點數（比如24個點，點越多越圓，最少3個點），返回圓的WKT（Well Known Text）文字
，WKT圖形繪製預覽工具：https://xiangyuecn.gitee.io/areacity-jsspider-statsgov/assets/geo-echarts.html
**/
static public String CreateSimpleCircleWKT(double lng, double lat, double radius, int pointCount) {
	double[][] points=CreateSimpleCircle(lng, lat, radius, pointCount);
	DecimalFormat df=new DecimalFormat("0.######");
	StringBuilder wkt=new StringBuilder("POLYGON((");
	for(int i=0;i<points.length;i++) {
		if(i>0)wkt.append(",");
		wkt.append(df.format(points[i][0])+" "+df.format(points[i][1]));
	}
	wkt.append("))");
	return wkt.toString();
}

js版原始碼

//測試：計算天壇到天安門的距離
console.log(Distance(116.410622, 39.881773, 116.397476, 39.908647));
//測試：生成天壇1公里範圍的圓形面
console.log(CreateSimpleCircleWKT(116.410622, 39.881773, 1000, 24));

/** 計算兩個座標的距離，單位米 **/
function Distance(lng1, lat1, lng2, lat2) {
    //採用Haversine formula演演算法，高德地圖的js計算程式碼，比較簡潔 https://www.cnblogs.com/ggz19/p/7551088.html
    var d=Math.PI/180;
    var f=lat1*d, h=lat2*d;
    var i=lng2*d - lng1*d;
    var e=(1 - Math.cos(h - f) + (1 - Math.cos(i)) * Math.cos(f) * Math.cos(h)) / 2;
    return 2 * 6378137 * Math.asin(Math.sqrt(e));
}

/** 以座標點為中心，簡單粗略的建立一個指定半徑的圓，半徑單位米，pointCount為構建圓的座標點數（比如24個點，點越多越圓，最少3個點），返回構成圓的座標點陣列 **/
function CreateSimpleCircle(lng, lat, radius, pointCount){
	//球面座標不會算，轉換成三角座標簡單點，經度代表值大約：0.01≈1km 0.1≈10km 1≈100km 10≈1000km
	var km=radius/1000;
	var a=km<5?0.01 :km<50?0.1 :km<500?1 :10;
	var b=Distance(lng, lat, lng+a, lat);
	var c=Distance(lng, lat, lng, lat+a);
	var rb=radius/b*a;
	var rc=radius/c*a;
	var arr=[];
	var n=0,step=360.0/pointCount,N=360-step/2; //注意浮點數±0.000000001的差異
	for(var i=0;n<N;i++,n+=step){
		var x=lng+rb*Math.cos(n*Math.PI/180);
		var y=lat+rc*Math.sin(n*Math.PI/180);
		arr[i]=[x, y];
	}
	arr.push([arr[0][0], arr[0][1]]); //閉環
	return arr;
};

/**
以座標點為中心，簡單粗略的建立一個指定半徑的圓，半徑單位米，pointCount為構建圓的座標點數（比如24個點，點越多越圓，最少3個點），返回圓的WKT（Well Known Text）文字
，WKT圖形繪製預覽工具：https://xiangyuecn.gitee.io/areacity-jsspider-statsgov/assets/geo-echarts.html
**/
function CreateSimpleCircleWKT(lng, lat, radius, pointCount){
	var points=CreateSimpleCircle(lng, lat, radius, pointCount);
    var wkt=["POLYGON(("];
    for(var i=0;i<points.length;i++) {
		wkt.push((i>0?",":"")+(+points[i][0].toFixed(6))+" "+(+points[i][1].toFixed(6)));
    }
    wkt.push("))");
    return wkt.join("");
};

線上繪製預覽效果

生成了圓形面的WKT文字後，可以貼上進線上預覽頁面繪製顯示：https://xiangyuecn.gitee.io/areacity-jsspider-statsgov/assets/geo-echarts.html ，方便程式碼偵錯，地圖上有測距功能，可以測量圓面的準確度；頁面上的畫圓功能，採用的就是js版的程式碼。

關於計算的精確度

兩個經緯度座標的距離計算，採用的Haversine formula演演算法，下面的計算程式碼中高德地圖的api同樣是Haversine formula演演算法，和百度地圖的計算結果誤差在0.2%以內：

bMap=window.BMapGL&&BMapGL.Map.prototype||{getDistance:function(a,b){return BMap.Map.prototype.getDistance(new BMap.Point(a.lng,a.lat),new BMap.Point(b.lng,b.lat)) }};

//地圖api計算【緯度】之間的距離，每一度之間的距離是相同的
bMap.getDistance({lng:111,lat:15},{lng:111,lat:16})        //百度111194.86米
new AMap.LngLat(111,15).distance(new AMap.LngLat(111,16)) //高德111319.49米

bMap.getDistance({lng:121,lat:55},{lng:121,lat:56})        //百度111194.78米
new AMap.LngLat(121,55).distance(new AMap.LngLat(121,56)) //高德111319.49米

//地圖api計算【經度】之間的距離，會隨著緯度的不同而不同
bMap.getDistance({lng:111,lat:15},{lng:112,lat:15})        //百度107405.91米
new AMap.LngLat(111,15).distance(new AMap.LngLat(112,15)) //高德107526.28米

bMap.getDistance({lng:111,lat:55},{lng:112,lat:55})        //百度63778.21米
new AMap.LngLat(111,55).distance(new AMap.LngLat(112,55)) //高德63849.69米
        //經度在相同緯度下每一度之間的距離是相同的
        bMap.getDistance({lng:121,lat:55},{lng:122,lat:55})        //百度63778.21米
        new AMap.LngLat(121,55).distance(new AMap.LngLat(122,55)) //高德63849.69米

對於構成圓的座標點的計算，使用的以前壓箱底的程式碼，計算比較簡單，經度和緯度分別計算出一度的距離長度，在等比例的換算出半徑對應的度數大小，比如算出來的經度1°是100km，那麼20km半徑對應的度數就是1° * 20 / 100 = 0.2°。

  • 緯度每一度之間的距離都是固定的長度，等比例換算後的結果是沒有誤差的。
  • 經度每一度之間的距離長度會隨著緯度變化，低緯度長高緯度短，圓面的上下緯度不同，半徑小的時候誤差小，半徑越大誤差越大。

所以按等比例換算在緯度上沒有問題，但經度上會產生一定的誤差，但只要圓的半徑不超過10km，誤差就能控制在0.5%以內，可通過下面程式碼直接計算觀察到：

//經度之間的距離，不同緯度下的誤差計算
d1=new AMap.LngLat(111,15.0).distance(new AMap.LngLat(111.2,15.0));
d2=new AMap.LngLat(111,15.2).distance(new AMap.LngLat(111.2,15.2));
console.log(d2, (d1-d2)/d2*100+"%") //低緯度下經度0.2°距離≈20km，緯度0.2°導致0.09%誤差
d1=new AMap.LngLat(111,55.0).distance(new AMap.LngLat(111.2,55.0));
d2=new AMap.LngLat(111,55.2).distance(new AMap.LngLat(111.2,55.2));
console.log(d2, (d1-d2)/d2*100+"%") //高緯度下經度0.2°距離≈10km，緯度0.2°導致0.50%誤差

d1=new AMap.LngLat(111,15).distance(new AMap.LngLat(112,15));
d2=new AMap.LngLat(111,16).distance(new AMap.LngLat(112,16));
console.log(d2, (d1-d2)/d2*100+"%") //低緯度下經度1°距離≈100km，緯度1°導致0.49%誤差
d1=new AMap.LngLat(111,55).distance(new AMap.LngLat(112,55));
d2=new AMap.LngLat(111,56).distance(new AMap.LngLat(112,56));
console.log(d2, (d1-d2)/d2*100+"%") //高緯度下經度1°距離≈60km，緯度1°導致2.57%誤差

由於是通過圓上的座標點連起來得到的一個圓，本質上是一個近似圓的多邊形，只要座標點足夠多就越接近圓，當座標點少的情況下肉眼可見的不是那麼圓（最少3個座標點，三角形），誤差也會很大。但考慮到點數越多，會導致使用上很多地方的計算量會變的很大，所以構成圓的座標點數也是一個綜合考慮的數量，我選擇24個座標點構成一個圓，每個象限6個點，點之間角度為15°。

【完】