根据图片中心点经纬度,计算图片四角的经纬度,并添加到地图中

 leaflet 详情总目录:传送 

        如何把图片在地图中显示?这里分两种情况:1、格式为tif,可以使用 L.imageOverlay,这个是使用geoserver的栅格数据展示图片,具体方法可以查看ImageMosaic 这篇文章;2、png/jpg等格式,如果知道图片对角经纬度,可以直接使用 L.imageOverlay 直接添加即可;如果仅有图片中心点经纬度,这里就需要通过计算来获取图片对角经纬度;如果没有图片相关的经纬度,那就放弃吧!

先看效果图:

今天主要是说一下根据图片中心点经纬度计算图片对角经纬度:

需要计算图片对角经纬度,就需要一些前提条件,如果没办法满足,这篇文章对你也就没啥作用了,不看也罢。

前提条件:1、相机传感器尺寸;2、相机拍摄分辨率;3、相机镜头焦距;4、拍摄高度。

计算公式:GSD = H * a / f;a = 传感器 / 分辨率; 

公式解释:H:拍摄高度 (m);f:相机镜头焦距 (mm);GSD:图片中每像素的真实距离 (m)

例如:传感器:35.7mm x 23.8mm ;拍摄分辨率:9504x6336;相机镜头焦距:40mm;拍摄高度:1000m;则:a = 35.7mm / 9504 = 0.00376mm(23.8mm / 6336 = 0.00376mm);GSD = 1000m * 0.00376mm / 40mm = 0.094m。(这些参数可以查看相机参数,或者问一下购买相机的客服)

前提条件已经准备好了,接下来直接上代码:

        /**
		 * 计算图片对角经纬度
		 * center 中心的坐标 
		 * latLength 图片高度像素
		 * lngLength 图片宽度像素
		 * gsd 
		 */
		computeLatlngByCenter(center, gsd, latLength, lngLength) {
			const latlng = center;
			let width = latLength * gsd;
			let height = lngLength * gsd;
			let one_lng_length = ((40076 * 1000) / 360) * Math.cos((latlng[0] * Math.PI) / 180);
            let one_lat_length = 111320 
			let lngOffSet = width / 2 / one_lng_length; 
			let latOffSet = height / one_lat_length / 2; 
			let north = latlng[0] + latOffSet;
			let south = latlng[0] - latOffSet;
			let east = latlng[1] + lngOffSet;
			let west = latlng[1] - lngOffSet;
			return [
				[north, west],
				[south, east],
			];
		},

注释:one_lng_length :1经度对应的纬线长度;one_lat_length :1纬度对应的经线长度,这个值是固定的;40076:赤道周长。以上数值都是在地球是圆形的基础上,和实际情况是有差距的,如果需要非常精确的数值,那这个方法就不合适了。

按照上面方法就可以算出图片对角的经纬度了,接下来就是把图片显示在地图上了:

        addImageLayer() {
			this.imageList = [
				{
					url: '/src/assets/image/A_01_DSC00520.jpg',
					center: [29.048664363492595, 119.76427712337504],
					height: 1000,
                    width:6336,
                    height:9504
				},
			];
			this.imageList.forEach((item) => {
				let gsd = 0.092;
				let imageBounds = this.computeLatlngByCenter(item.center, gsd,item.width,item.height);
				L.imageOverlay(item.url, imageBounds).addTo(this.map);
			});
		},

<think>好的,我现在需要解决用户的问题:如何根据中心点经纬度计算栅格的边界经纬度。首先,我得理解用户的需求。用户提到了参考引用中的内容,特别是引用[2]中提到的将经纬度转换为米坐标的方法,以及引用[1]中提到的计算几何的方法。看起来用户可能是在处理地理数据,比如栅格或矢量数据,需要根据中心点确定一个正方形的范围。 首先,我需要回忆经纬度到平面坐标的转换方法。通常,这种转换涉及到投影,比如Web墨卡托投影(EPSG:3857),这也是引用[2]中使用的方法。用户提到的lonLat2Mercator_函数就是将经纬度转换为米坐标的。这一步是关键,因为直接在地理坐标系(经纬度)下计算距离和范围会比较复杂,尤其是在高纬度地区,经度的距离会随着纬度变化而变化。 接下来,用户的问题是基于中心经纬度计算栅格的边框经纬度。假设栅格是一个正方形,边长为s米,那么需要找到四个边界的经纬度。步骤可能如下: 1. 将中心点经纬度转换为墨卡托坐标(米为单位)。 2. 根据边长s计算四个边界的米坐标。例如,中心点坐标是(x_center, y_center),那么左边界x_left = x_center - s/2,右边界x_right = x_center + s/2,同理上下边界的y坐标。 3. 将这四个边界的米坐标再转换回经纬度。 但需要注意,边长s是正方形的边长,转换为米后,需要确保在投影后的坐标系中,正方形的边确实对应实际地面上的等长。不过,由于墨卡托投影在高纬度地区会有变形,这种方法可能在远离赤道的地区不够准确,但对于小范围的栅格可能影响不大。 另外,用户提供的引用[2]中的代码似乎有一个问题:在计算y的时候,可能有一个错误。原代码中的第9行,计算y时使用了Math.tan的参数,正确的公式应该是将纬度转换为弧度,然后计算墨卡托投影的y坐标。正确的墨卡托投影公式应该是: y = ln(tan(π/4 + φ/2)),其中φ是纬度(以弧度为单位)。然后乘以地球半径(通常为6378137米),但在引用[2]中的代码中,地球周长被近似为20037508.34米(这是地球赤道的周长的一半,即赤道半径为6378137米,周长2πr约为40075016.68米,所以一半是20037508.34米)。所以,正确的转换应该是将经度直接乘以这个数除以180,而纬度则需要用上述公式转换。 不过用户提供的代码可能有些问题,例如变量y在第9行被计算后,第10行又乘以20037508.34/180,这可能正确,但需要验证。比如,假设输入的纬度是lat,转换为弧度后,计算tan((90 + lat) * π / 360)实际上等同于tan(π/4 + lat*π/360),因为90度对应π/2弧度,所以(90 + lat)度转换为弧度是(90 + lat)*π/180 / 2 = π/4 + lat*π/360。不过,正确的墨卡托投影公式应该是ln(tan(π/4 + φ/2)),其中φ是纬度弧度。所以如果lat是以度数表示的纬度,那么φ = lat * π/180,所以正确的表达式应为ln(tan(π/4 + (lat * π/180)/2)) = ln(tan(π/4 + lat * π/360)),这和代码中的第9行是一致的。因此,代码中的计算是正确的。然后第10行将y乘以20037508.34/180,相当于将结果缩放到米坐标,因为20037508.34米是赤道周长的一半,即每度的经度在赤道上的距离约为111319.49米(20037508.34/180),所以这可能适用于墨卡托投影的x和y坐标的转换。 接下来,计算边界的步骤需要详细展开: 假设中心点经纬度是(lon_center, lat_center),边长是s米。首先将中心点转换为墨卡托坐标(x_center, y_center)。然后,栅格在墨卡托坐标系中的范围是x_min = x_center - s/2,x_max = x_center + s/2,y_min = y_center - s/2,y_max = y_center + s/2。这样得到一个正方形的区域。然后将这四个点的坐标转换回经纬度,即得到对应的边界经纬度。 不过,这里需要注意的是,墨卡托投影是等的,但不是等面积的,且在高纬度地区y方向的拉伸会比较严重。因此,当s较大时,尤其是在高纬度地区,这种方法得到的“正方形”在实际地理空间中的形状可能不是真正的正方形,而是会有所变形。但对于小范围或低纬度地区,这种近似是可以接受的。 另外,用户可能在处理栅格数据时,需要根据中心点和像素大小(分辨率)来计算整个栅格的范围。例如,如果栅格的宽度是n像素,每个像素的分辨率是r米,那么边长s = n * r。然后按照上述方法计算边界。 现在,我需要将这些步骤转化为具体的公式和代码示例。同时,需要确保转换函数的正确性,参考用户提供的代码,但可能需要修正其中的潜在问题。 例如,用户提供的代码中的第9行: var y = Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180) 这里,lat是纬度,假设单位是度数。计算tan的参数是(90 + lat)度转换为弧度,即乘以π/360。但根据墨卡托投影的公式,正确的转换应该是: y = ln(tan(π/4 + φ/2)),其中φ是纬度弧度。φ = lat * π/180,所以π/4 + φ/2 = π/4 + (lat * π/180)/2 = π/4 + lat * π/360 = (90 + lat) * π/180 / 2 = (90 + lat) * π/360。所以代码中的这部分是正确的。然后,Math.log的结果是自然对数,之后除以(π/180)实际上将结果转换为度数?这可能有问题。因为墨卡托投影的y坐标应该是一个米单位的坐标。可能用户提供的代码存在错误,正确的墨卡托投影转换应该不需要除以(π/180),而是直接将结果乘以地球半径。 例如,正确的墨卡托投影的y坐标公式是: y = ln(tan(π/4 + φ/2)) * R 其中R是地球半径(约6378137米)。而用户代码中的处理是: y = ln(tan((90 + lat) * π / 360)) → 这等于ln(tan(π/4 + φ/2)),其中φ = lat * π/180。然后除以(π/180)得到的结果单位是度数?然后乘以20037508.34/180,这相当于将度数转换为米,因为20037508.34米是赤道半周长,即每度经度在赤道的长度约为20037508.34 * 2 / 360 ≈ 111319.49米。但这里可能混淆了经度和纬度的转换。正确的做法应该是经度的转换是直接乘以R * π/180,而纬度的转换如上述公式。 可能用户提供的代码中的转换存在错误,导致y坐标计算不准确。例如,正确的墨卡托投影的y坐标应该是: y = R * ln(tan(π/4 + φ/2)) 所以在代码中,应该直接计算这个值,而不需要除以(π/180)。但用户代码中的第9行和第10行: 第9行:y = Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180) 这实际上是将结果转换为度数?因为原计算得到的是以弧度为单位的数值,除以(π/180)会将其转换为度数。然后第10行乘以20037508.34/180,相当于将度数转换为米,乘以每度的米数。但这里可能存在问题,因为墨卡托投影的y坐标应该是直接以米为单位,不需要先转成度数再乘以每度的米数。正确的转换应该是: y = R * ln(tan(π/4 + φ/2)) 其中R=6378137米。所以对于纬度lat,φ = lat * π/180,则: y = 6378137 * ln(tan(π/4 + lat * π/360)) 然后,用户代码中的处理是: y = ln(tan((90 + lat)*π/360)) → 正确等同于ln(tan(π/4 + φ/2)),其中φ=lat*π/180 然后除以(π/180)得到的数值单位是弧度转换为度数?这里可能存在问题。例如,假设lat=0,则: y = ln(tan(π/4 + 0)) = ln(1) = 0 → 除以(π/180)后是0,乘以20037508.34/180后还是0,正确,因为赤道的y坐标是0。 但若lat=45度: (90+45)=135 → 135 * π/360 = 3π/8 → tan(3π/8) ≈ 2.4142 → ln(2.4142) ≈ 0.8814 → 除以π/180 ≈ 0.8814 / 0.01745 ≈ 50.5 → 然后乘以20037508.34/180 ≈ 50.5 * 111319.49 ≈ 5,620,000米。而根据正确的墨卡托公式,y=6378137 * ln(tan(π/4 + 45π/360)) = 6378137 * ln(tan(3π/8)) ≈ 6378137 * 0.8814 ≈ 5,622,000米,接近用户代码中的结果。所以用户提供的代码虽然在处理上有些绕,但结果可能是正确的,因为最终结果与正确公式得到的数值相近。因此,可能用户的代码是正确的,可以继续使用。 接下来,编写转换函数后,计算边界的步骤: 例如,中心点经纬度(lon_center, lat_center)转换为墨卡托坐标(x_center, y_center)。然后,每个边界的墨卡托坐标为: x_min = x_center - s/2 x_max = x_center + s/2 y_min = y_center - s/2 y_max = y_center + s/2 然后将这四个点转换回经纬度,得到四个点的经纬度,但实际只需要最大最小的经纬度即可。例如,左边界对应x_min,右边界x_max,上下边界y_min和y_max对应的纬度。不过,需要注意的是,在墨卡托投影中,x对应经度,y对应纬度,所以将x_min和x_max转换回经度时,直接使用反函数,而y_min和y_max转换回纬度时,需要使用反函数。 具体来说,反转换函数需要将墨卡托坐标的x和y转换为经度和纬度: 对于x坐标,经度lon = (x * 180) / 20037508.34 对于y坐标,需要解方程: y_merc = 20037508.34 / 180 * [ ln(tan( (90 + lat) * π / 360 )) ] 所以反函数为: lat = ( arctan( exp( y_merc * 180 / 20037508.34 * π / 180 ) ) - π/4 ) * 2 * 180 / π 或者简化步骤: 从y_merc得到: y_merc = (20037508.34 / 180) * [ ln(tan( (90 + lat) * π / 360 )) ] 因此: ln(tan( (90 + lat) * π / 360 )) = y_merc * 180 / 20037508.34 取指数: tan( (90 + lat) * π / 360 ) = exp( y_merc * 180 / 20037508.34 ) 然后求反正切: (90 + lat) * π / 360 = arctan( exp( y_merc * 180 / 20037508.34 ) ) 解得: lat = [ 360 * arctan( exp( y_merc * 180 / 20037508.34 ) ) / π - 90 ] 或者化简: lat = ( 2 * arctan( exp( y_merc * π / 20037508.34 ) ) - π/2 ) * 180 / π 这可能更简洁。例如,对于y_merc的米坐标,转换为纬度: lat = (2 * arctan(exp( y_merc * π / 20037508.34 )) - π/2) * 180 / π 现在,将这些步骤整理成公式,写出对应的代码示例。 总结步骤: 1. 定义转换函数,将经纬度转换为墨卡托坐标(米)。 2. 计算中心点的墨卡托坐标。 3. 根据边长计算边界的墨卡托坐标4. 将边界的墨卡托坐标转换回经纬度。 需要注意的是,用户可能需要处理的是正方形的边长为s米,那么在墨卡托坐标系中,正方形的边长在x和y方向都是s米。但由于墨卡托投影的特性,这在实际地理中可能非精确的正方形,尤其是在高纬度地区。 现在,编写具体的公式和代码示例: 转换函数参考用户提供的代码,调整可能的错误: 例如,修复用户代码中的var y在第9行被重新赋值的问题,可能需要使用val: def lonLat2Mercator(lon: Double, lat: Double): (Double, Double) = { val x = lon * 20037508.34 / 180 val y = Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180) val y_merc = y * 20037508.34 / 180 (x, y_merc) } 但根据之前的分析,用户的代码可能已经正确,因为第9行得到的是一个以度数为单位的值,然后乘以20037508.34/180转换为米。或者,可能更直接的方式是: 正确的墨卡托投影的转换函数应该是: def lonLat2Mercator(lon: Double, lat: Double): (Double, Double) = { val x = lon * Math.PI / 180 * 6378137 val y = Math.log( Math.tan( Math.PI/4 + Math.toRadians(lat)/2 )) * 6378137 (x, y) } 但这样计算出来的x和y会以米为单位。然而,用户提供的代码中使用了20037508.34这个值,这实际上是赤道周长的一半,即赤道半径为6378137米,赤道周长为2 * π * 6378137 ≈ 40075016.68米,所以20037508.34是半周长。因此,经度每度的长度在赤道是20037508.34 * 2 / 360 ≈ 111319.49米。用户代码中将经度lon乘以20037508.34/180,等同于lon * 111319.49,这与直接使用米为单位的墨卡托投影一致。因此,用户提供的代码是正确的,可能只是另一种表达方式。 因此,转换函数可以保留用户的实现,然后反函数则需要根据同样的参数来推导。 例如,反函数mercator2LonLat: def mercator2LonLat(x: Double, y: Double): (Double, Double) = { val lon = x * 180 / 20037508.34 val temp = y * 180 / 20037508.34 val lat = Math.atan( Math.exp( temp * Math.PI / 180 )) * 360 / Math.PI - 90 (lon, lat) } 这里,temp是将y坐标转换为度数,因为用户的正向转换中,y坐标是: y_merc = [ ln(tan(...)) / (π/180) ] * (20037508.34/180 ) 即 y_merc = [ln(...) * 180 / π ] * (20037508.34 / 180 ) = ln(...) * 20037508.34 / π 但可能反函数的推导需要更精确。例如,假设正向转换是: y_merc = ( ln(tan( (90 + lat) * π/360 )) ) * (20037508.34 / 180 ) 那么,反转换时: y_merc * 180 / 20037508.34 = ln(tan( (90 + lat) * π/360 )) 取指数两边: exp( y_merc * 180 / 20037508.34 ) = tan( (90 + lat) * π/360 ) 然后取反正切: arctan( exp( y_merc * 180 / 20037508.34 ) ) = (90 + lat) * π/360 解出lat: lat = ( arctan( exp( y_merc * 180 / 20037508.34 )) * 360 / π ) - 90 所以反函数中的纬度计算应为: val lat = Math.atan( Math.exp( y * 180 / 20037508.34 )) * 360 / Math.PI - 90 而经度lon的计算直接是: lon = x * 180 / 20037508.34 这样,反函数就正确了。 现在,步骤总结: 1. 将中心点(lon_center, lat_center)转换为墨卡托坐标(x_center, y_center)。 2. 计算边界的墨卡托坐标:x_min = x_center - s/2,x_max = x_center + s/2,y_min = y_center - s/2,y_max = y_center + s/2。 3. 将这四个点(x_min, y_min), (x_max, y_max)等转换回经纬度,得到对应的边界经纬度。 不过,实际需要的可能是最大最小的经度和纬度,所以需要将x_min和x_max转换为经度的左右边界,y_min和y_max转换为纬度的上下边界。 例如,左边界经度是x_min对应的lon_left = x_min * 180 / 20037508.34,右边界lon_right = x_max * 180 / 20037508.34。同样,下边界纬度lat_bottom = mercator2Lat(y_min),上边界lat_top = mercator2Lat(y_max)。 因此,最终的边界经纬度为: 最小经度:lon_left = x_center - s/2 转换后的经度 最大经度:lon_right = x_center + s/2 转换后的经度 最小纬度:lat_bottom = y_center - s/2 转换后的纬度 最大纬度:lat_top = y_center + s/2 转换后的纬度 但需要注意,在墨卡托投影中,y轴的方向可能与纬度增长方向一致,因此y_center + s/2对应北边的纬度更高,而y_center - s/2对应南边的纬度更低。 现在,将这些步骤用具体的例子说明。假设中心点经度是116.0,纬度是39.0,边长s是1000米。那么: 1. 转换中心点到墨卡托坐标: 使用用户提供的函数: x_center = 116.0 * 20037508.34 / 180 ≈ 116 * 111319.49 ≈ 12,913, 000 米(具体数值需要计算) y_center的计算: lat = 39.0,代入函数: y = ln(tan((90 + 39) * π/360 )) / (π/180) → (129 * π/360) ≈ 1.117弧度,tan≈ 2.0 → ln(2.0) ≈ 0.6931 → 0.6931 / (π/180) ≈ 0.6931 * 57.2958 ≈ 39.7 → y_merc = 39.7 * 20037508.34 / 180 ≈ 39.7 * 111319.49 ≈ 4,421, 000 米 2. 计算x_min = x_center - 500米,x_max = x_center + 500米;y_min = y_center -500米,y_max = y_center +500米. 3. 转换这些点回经纬度: 例如,x_min的经度: (x_center -500) * 180 / 20037508.34 这将得到左右边界的经度,上下边界的纬度同理。 最后,用户可能需要将这些步骤写成代码,例如用Python实现,或者使用其他语言。 不过用户提供的引用[3]提到geopandas可以方便处理地理数据,可能用户更倾向于使用现有的库函数,例如使用pyproj库进行坐标转换,而不是手动实现。例如,使用pyproj的Transformer来转换坐标系。 例如,使用EPSG:4326(WGS84经纬度)和EPSG:3857(Web墨卡托)之间的转换: from pyproj import Transformer transformer_to_mercator = Transformer.from_crs("EPSG:4326", "EPSG:3857") transformer_to_lonlat = Transformer.from_crs("EPSG:3857", "EPSG:4326") # 中心点经纬度 lon_center = 116.0 lat_center = 39.0 # 转换为墨卡托坐标 x_center, y_center = transformer_to_mercator.transform(lat_center, lon_center) # 计算边界的墨卡托坐标 s = 1000 # 边长1000米 x_min = x_center - s/2 x_max = x_center + s/2 y_min = y_center - s/2 y_max = y_center + s/2 # 转换墨卡托坐标经纬度 lat_min, lon_min = transformer_to_lonlat.transform(y_min, x_min) lat_max, lon_max = transformer_to_lonlat.transform(y_max, x_max) # 或者更准确地说,转换四个点: # 左下 (x_min, y_min) → (lat_min, lon_min) # 右上 (x_max, y_max) → (lat_max, lon_max) # 但需要注意,墨卡托坐标的x和y对应的是经度和纬度,但转换时可能需要调整参数的顺序(通常纬度在前,经度在后) 不过,pyproj的transform方法参数顺序通常是纬度、经度,或者根据坐标系的定义。例如,在EPSG:4326中,坐标顺序是纬度、经度,而EPSG:3857是东向、北向。因此,转换时需要注意参数的顺序。 例如,正确的转换可能是: x_center, y_center = transformer_to_mercator.transform(lat_center, lon_center) → 因为EPSG:4326的参数顺序是纬度、经度,而EPSG:3857是东向(经度对应的x)、北向(纬度对应的y)。 因此,正确的转换应该是: x_center, y_center = transformer_to_mercator.transform(lon_center, lat_center) → 这里可能我之前的参数顺序有误,需要根据pyproj的文档确认。通常,对于EPSG:4326,参数顺序是经度、纬度,但有些情况下可能要求纬度、经度。需要仔细检查。 例如,根据pyproj的文档,默认的EPSG:4326使用的是纬度、经度顺序,但可以通过设置always_xy=True来使用经度、纬度顺序。所以,更可靠的代码是: transformer_to_mercator = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) x_center, y_center = transformer_to_mercator.transform(lon_center, lat_center) 这样,transform方法的参数顺序就是经度、纬度。 因此,正确的转换步骤应该是: 1. 使用pyproj的Transformer,设置always_xy=True以确保参数顺序为经度、纬度。 2. 转换中心点到墨卡托坐标。 3. 计算边界的墨卡托坐标4. 转换回经纬度,得到边界的经度和纬度。 这样,可以避免手动实现转换函数中的潜在错误,特别是墨卡托投影的反函数较为复杂。 因此,最终的解决方案可能建议用户使用现有的库如pyproj来进行坐标转换,这样更准确和简便,同时引用用户提供的参考中的方法。 总结来说,用户的解决方案有两种途径: 1. 手动实现经纬度到墨卡托坐标的转换及反转换,根据用户提供的代码进行调整。 2. 使用现有的地理处理库(如pyproj、geopandas)来简化转换过程。 考虑到用户提到了引用[3]中的geopandas,可能用户更倾向于使用这些库。因此,在回答中应该同时介绍两种方法,手动公式和库函数的使用。 现在,整理这些思考过程,形成结构清晰的回答,包括步骤、公式、代码示例,注意引用用户提供的参考内容。</think>根据中心点经纬度计算栅格边界的经纬度,通常需要借助地理坐标系与投影坐标系的转换。以下是具体方法和公式说明: --- ### **步骤说明** 1. **将中心点经纬度转换为墨卡托投影坐标(米)** 使用Web墨卡托投影(EPSG:3857),将经纬度转换为平面坐标,便于进行米制单位的距离计算[^2]。 转换公式: $$ x = \frac{lon \cdot 20037508.34}{180} $$ $$ y = \frac{\ln\left(\tan\left(\frac{(90 + lat) \cdot \pi}{360}\right)\right) \cdot 20037508.34}{\pi \cdot 180} $$ 2. **计算栅格边界的墨卡托坐标** 假设栅格边长为$s$米,中心点墨卡托坐标为$(x_c, y_c)$,则边界坐标为: $$ x_{\text{min}} = x_c - \frac{s}{2}, \quad x_{\text{max}} = x_c + \frac{s}{2} $$ $$ y_{\text{min}} = y_c - \frac{s}{2}, \quad y_{\text{max}} = y_c + \frac{s}{2} $$ 3. **将边界墨卡托坐标转换回经纬度** 反向转换公式: $$ lon_{\text{boundary}} = \frac{x_{\text{boundary}} \cdot 180}{20037508.34} $$ $$ lat_{\text{boundary}} = \frac{2 \cdot \arctan\left(\exp\left(\frac{y_{\text{boundary}} \cdot \pi}{20037508.34}\right)\right) - \frac{\pi}{2}} \cdot 180}{\pi} $$ --- ### **代码实现(Python)** #### 方法1:手动计算 ```python import math def lonlat_to_mercator(lon, lat): x = lon * 20037508.34 / 180 y = math.log(math.tan((90 + lat) * math.pi / 360)) * 20037508.34 / (math.pi * 180) return x, y def mercator_to_lonlat(x, y): lon = x * 180 / 20037508.34 lat = (2 * math.atan(math.exp(y * math.pi / 20037508.34)) - math.pi/2) * 180 / math.pi return lon, lat # 示例:中心点经纬度 (116.0, 39.0),边长 1000 米 lon_center, lat_center = 116.0, 39.0 s = 1000 # 转换为墨卡托坐标 x_center, y_center = lonlat_to_mercator(lon_center, lat_center) # 计算边界 x_min = x_center - s/2 x_max = x_center + s/2 y_min = y_center - s/2 y_max = y_center + s/2 # 转换回经纬度 lon_min, lat_min = mercator_to_lonlat(x_min, y_min) lon_max, lat_max = mercator_to_lonlat(x_max, y_max) ``` #### 方法2:使用`pyproj`库 ```python from pyproj import Transformer # 定义坐标系转换器 transformer_to_mercator = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) transformer_to_lonlat = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True) # 计算中心点墨卡托坐标 x_center, y_center = transformer_to_mercator.transform(lon_center, lat_center) # 计算边界转换回经纬度 lon_min, lat_min = transformer_to_lonlat.transform(x_min, y_min) lon_max, lat_max = transformer_to_lonlat.transform(x_max, y_max) ``` --- ### **注意事项** 1. **投影变形**:墨卡托投影在高纬度地区会导致较大的形变,建议在小范围(如城市级)使用此方法[^2]。 2. **库函数推荐**:使用`geopandas`或`pyproj`可简化投影转换流程[^3]。 ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

travel_wsy

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值