python - GeoTIFF 数据文件中的纬度/经度(公里)与计算的距离不匹配
问题描述
我正在使用gdalinfo
基于THIS从 GeoTIFF 文件中提取信息,我得到以下信息:
Driver: GTiff/GeoTIFF
Files: myfile.tif
myfile.tif.aux.xml
Size is 953, 2824
Coordinate System is:
PROJCS["Lambert_Conformal_Conic_2SP_World_Geodetic_System_1984",
GEOGCS["GCS_World_Geodetic_System_1984",
DATUM["D_World_Geodetic_System_1984",
SPHEROID["WGS_1984",6378137,298.257223563,
AUTHORITY["EPSG","7030"]]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",66],
PARAMETER["standard_parallel_2",78],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-42],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-432678.541899999952875,9726108.134099999442697)
Pixel Size = (300.000000000000000,-300.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -432678.542, 9726108.134) ( 53d58'12.29"W, 70d52'36.63"N)
Lower Left ( -432678.542, 8878908.134) ( 50d38' 9.28"W, 63d22'47.25"N)
Upper Right ( -146778.542, 9726108.134) ( 46d 6'31.44"W, 71d13'15.48"N)
Lower Right ( -146778.542, 8878908.134) ( 44d56'51.10"W, 63d37'31.84"N)
Center ( -289728.542, 9302508.134) ( 48d45'16.75"W, 67d18'24.75"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
Min=0.900 Max=1.100
Minimum=0.900, Maximum=1.100, Mean=0.993, StdDev=0.025
NoData Value=-3.4028234663852886e+38
Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=1.1000000238419
STATISTICS_MEAN=0.99293445610505
STATISTICS_MINIMUM=0.89999997615814
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=0.025099979310242
我正在尝试根据纬度/经度参考重现 4 个角的公里距离。我可以得到经度,但不知何故我的纬度偏离了……这是我的做法:
import geopy
from geopy import distance
from geopy.location import Point
p1 = Point('''70 52'36.63" N 53 58'12.29" W''')
p2 = Point('''0 N 53 58'12.29" W''') #Same longitude, but latitude set to 0
distance.distance(p1,p2).m # WGS84 distance in m
输出基础gdalinfo
,我希望能回来9726108.134
,但相反,我出去了7866807.760742959
。我认为这是一个投影错误,太过分了。关于我做错了什么的任何其他建议?
请注意,当我在此示例中复制 km 时,我可以使用我的方法正确获得以 km 为单位的纬度...
解决方案
您的距离检查做了两个假设:
- 点 0° N, 0° E 也在重新投影中的 0, 0 处。
- x 和 y 轴在两个投影中具有相同的对齐方式
然而,这两种说法都不成立。这是一个图表来说明两个预测之间的差异:
左图显示了在您的 GeoTiff 中描绘为 WGS84 中的绿色多边形的区域,右图在重新投影中显示了相同的内容。如您所见,点 0° N、0° E 现在位于 7633606.089886177、2779916.995372205(而重新投影的点 0、0 在 WGS84 中将位于 0° N、-42° E)。此外,在右侧形成完美矩形的图像区域在 WGS84 中被扭曲。
长话短说:由于预测固有的差异,您执行的检查不会得出您的预期结果。作为近似值,您可以改为使用四个角的坐标进行检查,结果将证明(大致)正确。这是使用左上角和左下角的示例:
p1 = geopy_Point('''70 52'36.63" N 53 58'12.29" W''')
p2 = geopy_Point('''63 22'47.25" N 50 38' 9.28" W''')
distance.distance(p1,p2).m
返回:
848195.724897564
这非常接近预期距离 (9726108.134 - 8878908.134 = 847200),并且由于给定 WGS84 坐标的四舍五入而略有偏差。
希望这可以帮助!
编辑:
以下是如何从 WGS84 转换为重新投影的建议:
import osr
from pyproj import Proj, transform
def reproject(x_in, y_in):
proj_in = Proj(init='epsg:4326')
osr_spat_ref = osr.SpatialReference()
osr_spat_ref.ImportFromWkt('''PROJCS["Lambert_Conformal_Conic_2SP_World_Geodetic_System_1984",
GEOGCS["GCS_World_Geodetic_System_1984",
DATUM["D_World_Geodetic_System_1984",
SPHEROID["WGS_1984",6378137,298.257223563,
AUTHORITY["EPSG","7030"]]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",66],
PARAMETER["standard_parallel_2",78],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-42],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]''')
spat_ref_pyproj = osr_spat_ref.ExportToProj4()
proj_out = Proj(spat_ref_pyproj)
out = transform(proj_in, proj_out, x_in, y_in)
return out
只需输入经度和纬度坐标作为x_in
和y_in
分别获得重新投影的坐标。从这些中,您可以通过扣除角点的坐标来确定您感兴趣的图像的哪个图块。
推荐阅读
- python - 执行 Sql 脚本返回“无”
- python - 为什么对象在 for 循环中加速?
- linux - 带有单独结果的 bash 简单搜索脚本
- r - 我想提高我的 hts 时间序列预测的准确性(MAPE)
- php - 递归循环遍历多维数组并对每个键值做一些事情
- jenkins-pipeline - 来自 AWS ECR 的 Jenkins 管道 Docker 代理
- c# - c# dllimport 试图加载格式不正确的程序
- c++ - set_union() 不适用于一组字符串
- django - 如何在 django 自定义保存方法中访问数据库对象?
- docusignapi - 如果文档是使用 Docusign API 签名的,如何访问已完成的文档?