假设有一张含有地理坐标信息的tif格式的图片及其对应的jpg或者png格式的普通图片
如下图所示：
其中第一张为tif格式的地理信息图，第二张为按照一定比例下采样之后转换得到的普通jpg图片
如何计算jpg图片中任意两点之间的实际距离呢？
比如，计算图中间位置的桥的长度，也就是桥两个端点之间的距离
在tif图片中，如果我们知道了点的经纬度就可以使用geopy包中的geodesic函数计算两点间的距离
代码如下
from geopy.distance import geodesic
distance = geodesic((lat1,lon1), (lat2,lon2)).meters
其中lat1，lon1为第一个点的纬度和经度
所以问题的关键在于如何找到jpg图片中的点对应的经纬度
为了达到该目标需要经过3个步骤
步骤一
找到jpg中的点在图片中所在的相对图片宽和高的比例
代码如下：
import cv2 as cv
image = cv.imread('example.jpg')
h, w, c = image.shape
# 选取两个点
pick1_x, pick1_y = 2022, 2160
pick2_x, pick2_y = 1902, 2160
pick1_x_ratio, pick1_y_ratio = pick1_x/w, pick1_y/h
pick2_x_ratio, pick2_y_ratio = pick2_x/w, pick2_y/h
步骤二
按照比例找到两个点对应的tif图中的行与列
代码如下：
import rasterio
with rasterio.open(tif_file) as src:
tmp_height = src.height # 数据高度
tmp_width = src.width # 数据宽度
pick1_col, pick1_row = int(pick1_x_ratio*tmp_width), int(pick1_y_ratio*tmp_height)
pick2_col, pick2_row = int(pick2_x_ratio*tmp_width), int(pick2_y_ratio*tmp_height)
步骤三
将行列号转换为以地图单位（Map units）表示的坐标
代码如下：
with rasterio.open(tif_file) as src:
pick1_lon_map, pick1_lat_map = src.xy(pick1_row, pick1_col)
pick2_lon_map, pick2_lat_map = src.xy(pick2_row, pick2_col)
步骤四
将地图单位坐标转换为十进制度数（Decimal Degrees）坐标，
此处应注意使用的坐标系
代码如下：
import pyproj
import rasterio
with rasterio.open(tif_file) as src:
# 定义源坐标系统和目标坐标系统（这里以 EPSG:4326 作为目标坐标系统，即 WGS84 经纬度坐标系统）
src_proj = pyproj.Proj(src.crs) # src.crs为获取到的tif图中的源坐标系统
dst_proj = pyproj.Proj(init='EPSG:4326')
print("picked points 坐标转换")
pick1_long, pick1_latt = pyproj.transform(src_proj, dst_proj, pick1_lon_map, pick1_lat_map)
pick2_long, pick2_latt = pyproj.transform(src_proj, dst_proj, pick2_lon_map, pick2_lat_map)
步骤五
根据十进制坐标，计算两点之间的距离
代码如下：
from geopy.distance import geodesic
# 纬度在前，经度在后
distance = geodesic((pick1_latt, pick1_long), (pick2_latt, pick2_long)).meters
print(f"选取的两点之间的距离：{distance} meters")
选取的两点为桥的两个端点，因此最终计算得到桥长为15.57米
整体代码如下：
import cv2 as cv
import rasterio
from geopy.distance import geodesic
import pyproj
# 打开图像文件
tif_file_path = 'path_to_tif_file'
jpg_file_path = 'path_to_jpg_file'
image = cv.imread(jpg_file_path)
h, w, c = image.shape
# 打开TIFF文件
with rasterio.open(tif_file_path) as src:
# 选取桥头两个点
pick1_x, pick1_y = 2022, 2160
pick2_x, pick2_y = 1902, 2160
tmp_height = src.height # 数据高度，行数
tmp_width = src.width # 数据宽度，列数
# step1 and step2: 计算两点所占比例及其在tif中对应的行列号
pick1_col, pick1_row = int(pick1_x*tmp_width/w), int(pick1_y*tmp_height/h)
pick2_col, pick2_row = int(pick2_x*tmp_width/w), int(pick2_y*tmp_height/h)
# step3: 计算对应的地图单位坐标
print("picked point 经纬度")
pick1_lon, pick1_lat = src.xy(pick1_row, pick1_col)
pick2_lon, pick2_lat = src.xy(pick2_row, pick2_col)
print("pick1经纬度:", pick1_lon, ' ', pick1_lat)
print("pick2经纬度:", pick2_lon, ' ', pick2_lat)
# 定义源坐标与目标坐标系统（这里以 EPSG:4326 作为目标坐标系统，即 WGS84 经纬度坐标系统）
src_proj = pyproj.Proj(src.crs)
dst_proj = pyproj.Proj(init='EPSG:4326')
# step 4: 将地图单位坐标转为十进制经纬度坐标
print("picked points 坐标转换")
pick1_long, pick1_latt = pyproj.transform(src_proj, dst_proj, pick1_lon, pick1_lat)
pick2_long, pick2_latt = pyproj.transform(src_proj, dst_proj, pick2_lon, pick2_lat)
print(f"第一个点经纬度：{pick1_long} {pick1_latt}")
print(f"第二个点经纬度：{pick2_long} {pick2_latt}")
# step5: 计算两点之间的距离
distance = geodesic((pick1_latt, pick1_long), (pick2_latt, pick2_long)).meters
print(f"选取的两点之间的距离：{distance} meters")
还有一个转换值得关注，
即行列号与地图单位坐标之间的互相转换
代码如下：
import rasterio
with rasterio.open(tif_file) as src:
# 行列号转地图单位坐标
lon_map, lat_map = src.xy(row, col)
# 地图单位坐标转行列号
row, col = src.index(lon_map, lat_map)
