Lat/lon GeoTIFF 元数据文件中以公里为单位的计算距离不匹配

Question

我基于 THIS 使用 gdalinfo 从 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

我正在尝试根据 lat/lon 度参考重现 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。我认为这是一个投影错误太多了。关于我做错了什么还有其他建议吗？

请注意，当我在 THIS EXAMPLE 中重现公里时，我可以使用我的方法正确获得以公里为单位的纬度...

Answer 1

您的距离检查有两个假设：

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

哪个returns:

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即可得到重投影坐标。从这些中，您可以通过扣除角的坐标来找出您感兴趣的图像的哪个图块。

Lat/lon GeoTIFF 元数据文件中以公里为单位的计算距离不匹配

Lat/lon in km in GeoTIFF medata file not matching calculated distance

python

tiff

geolocation

gdal

geotiff