摘要:依据经纬坐标计算两点距离并求图像分辨率输入图像文件输出两点距离和图像分辨率缺点未将投影方式纳入考虑范围,处理高纬地区可能精度失真。地球平均半径,用公式计算球面两点间的距离。
依据经纬坐标计算两点距离并求图像分辨率
输入:geotiff图像文件
输出:两点距离和图像分辨率
缺点:未将投影方式纳入考虑范围, 处理高纬地区可能精度失真。
import math, gdalfrom geopy.distance import geodesicfrom osgeo import gdal, osrimport numpy as npfrom math import sin, asin, cos, radians, fabs, sqrt EARTH_RADIUS = 6371 # 地球平均半径,6371kmdef hav(theta): s = sin(theta / 2) return s * s def get_distance_hav(lat0, lng0, lat1, lng1): """用haversine公式计算球面两点间的距离。""" # 经纬度转换成弧度 lat0 = radians(lat0) lat1 = radians(lat1) lng0 = radians(lng0) lng1 = radians(lng1) dlng = fabs(lng0 - lng1) dlat = fabs(lat0 - lat1) h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng) distance = 2 * EARTH_RADIUS * asin(sqrt(h)) return distancedef radi(d): return d * 3.1415926 / 180def read_img(filename): dataset = gdal.Open(filename) #打开文件 im_width = dataset.RasterXSize #栅格矩阵的列数 im_height = dataset.RasterYSize #栅格矩阵的行数 im_bands = dataset.RasterCount #波段数 im_geotrans = dataset.GetGeoTransform() #仿射矩阵,左上角像素的大地坐标和像素分辨率 im_proj = dataset.GetProjection() #地图投影信息,字符串表示 im_data = dataset.ReadAsArray(0,0,im_width,im_height) # del dataset print(im_bands, im_height, im_width, im_geotrans, im_proj) return im_bands, im_data, datasetdef getDist(lat1, lat2, lon1, lon2): radlat = radi(lat1) - radi(lat2) radlon = radi(lon1) - radi(lon2) print(radlat, radlon) dist = 2*math.asin(math.sqrt(math.pow(math.sin(radlat/2), 2) + math.cos(radi(lat1)) * math.cos(radi(lat2))*math.pow(math.sin(radlon/2), 2))) #math.pow dist = math.floor(dist * 6378137 * 10000) / 10000 return distdef getDist2(lat1, lat2, lon1, lon2): # dist = geodesic((30.28708, 120.12802999999997), (28.7427, 115.86572000000001)).km # dist = geodesic((lat1, lon1), (lat2, lon2)) return dist def getSRSPair(dataset): """ 获得给定数据的投影参考系和地理参考系 :param dataset: GDAL地理数据 :return: 投影参考系和地理参考系 """ prosrs = osr.SpatialReference() prosrs.ImportFromWkt(dataset.GetProjection()) geosrs = prosrs.CloneGeogCS() return prosrs, geosrsdef geo2lonlat(dataset, x, y): """ 将投影坐标转为经纬度坐标(具体的投影坐标系由给定数据确定) :param dataset: GDAL地理数据 :param x: 投影坐标x :param y: 投影坐标y :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat) """ prosrs, geosrs = getSRSPair(dataset) ct = osr.CoordinateTransformation(prosrs, geosrs) coords = ct.TransformPoint(x, y) return coords[:2]def resolu(filename): im_bands, im_data, dataset = read_img(filename) trans = dataset.GetGeoTransform() x0, y0, xResolution, yResolution = trans[0], trans[3], trans[1], trans[5] coords = geo2lonlat(dataset, x0, y0) print("(%s, %s)->(%s, %s)" % (trans[0], trans[3], coords[0], coords[1])) # y0, x0 = coords[0], coords[1] y1 = y0 - yResolution * im_data.shape[1] y0 = coords[0] x0 = coords[1] coords = geo2lonlat(dataset, x0, y1) print("(%s, %s)->(%s, %s)" % (x0, y1, coords[0], coords[1])) print(x0, y0, y1, xResolution, yResolution, im_data.shape) y1 = coords[0] print(y0, y1, x0, x0) ylength1 = getDist(y0, y1, x0, x0) ylength = getDist2(y0, y1, x0, x0) print("ylength, ylength1 :: ", ylength, ylength1) yResolu = ylength1 / im_data.shape[1] x1 = x0 + xResolution * im_data.shape[2] # xlength = getDist(y0, y0, x0, x1) xlength = getDist2(y0, y0, x0, x1) xResolu = xlength / im_data.shape[2] return xResolu, yResolu, xlength, ylengthif __name__ == "__main__": # 依据经纬坐标计算两点距离并求图像分辨率 xResolu, yResolu, xlength, ylength = resolu("D:/111/Rectangle_03_卫图/Rectangle_03_卫图_Level_18.tif") print("xResolu, yResolu :: ", xResolu, yResolu, xlength, ylength)
文章版权归作者所有,未经允许请勿转载,若此文章存在违规行为,您可以联系管理员删除。
转载请注明本文地址:https://www.ucloud.cn/yun/124472.html
摘要:最近工作需要,网上搜索了下根据经纬度计算两地距离的方法,发现要么是几何法,画图作一堆辅助线,然后证明推理,要么二话不说直接套公式。球体上两地的最短距离就是经过两点的大圆的劣弧长度。 最近工作需要,网上搜索了下根据经纬度计算两地距离的方法,发现要么是几何法,画图、作一堆辅助线,然后证明推理,要么二话不说直接套公式。这篇文章介绍一种容易理解的方式来求这个距离。 0b00 思路 地球是个不规...
摘要:霍夫变换霍夫变换是图像处理中的一种特征提取技术,它通过一种投票算法检测具有特定形状的物体。该过程在一个参数空间中通过计算累计结果的局部最大值得到一个符合该特定形状的集合作为霍夫变换的结果。 参考的一些文章以及论文我都会给大家分享出来 —— 链接就贴在原文,论文我上传到资源中去,大家可以免...
摘要:前言最近在帮朋友商家写小程序,所以看了不少关于小程序的知识,总结一下计算距离这条线。 前言 最近在帮朋友(商家)写小程序,所以看了不少关于小程序的知识,总结一下计算距离这条线。 思路 一共有两种方法,各有利弊:1.利用小程序的wx.getLocation 方法得到用户的经纬度,然后用已知的商家的经纬进行计算;2.利用腾讯地图位置服务calculateDistance直接计算 先熟悉下两...
摘要:计算精度与谷歌地图的距离精度差不多,相差范围在米以下。以上代码大部分来自网上收集,经过验证过的,可以使用 根据经纬度计算距离公式 showImg(https://segmentfault.com/img/bV6zX2?w=437&h=76); 图片来自互联网 对上面的公式解释如下: Lung1 Lat1表示A点经纬度, Lung2 Lat2表示B点经纬度; a=Lat1 – Lat2...
摘要:计算精度与谷歌地图的距离精度差不多,相差范围在米以下。以上代码大部分来自网上收集,经过验证过的,可以使用 根据经纬度计算距离公式 showImg(https://segmentfault.com/img/bV6zX2?w=437&h=76); 图片来自互联网 对上面的公式解释如下: Lung1 Lat1表示A点经纬度, Lung2 Lat2表示B点经纬度; a=Lat1 – Lat2...
阅读 3413·2021-11-24 10:22
阅读 3123·2021-11-23 09:51
阅读 3464·2021-11-22 09:34
阅读 3230·2021-11-19 09:40
阅读 2277·2021-11-15 11:39
阅读 1331·2021-10-14 09:42
阅读 3502·2021-10-08 10:04
阅读 1393·2019-08-30 15:52