python – 为什么Basemap南极立体地图投影坐标与同一投影中的数据集不一致？

2019年7月20日 498次阅读

一些基于卫星的地球观测产品提供纬度/经度信息,而其他产品提供给定网格投影内的X / Y坐标(并且还有一些具有两者,参见示例).

我在第二种情况下的方法是设置一个底图贴图,该贴图具有与数据提供者给定的相同参数(投影,椭圆体,贴图原点),其方式是给定的X / Y值等于底图坐标.但是,如果我这样做,地理定位与包括Basemap海岸线在内的其他数据集不一致.

我从不同的可靠来源获得了三种不同的数据集.对于最小的例子,我使用
U.S. Geological Survey提供的Landsat数据,其中包括南极立体网格的X / Y坐标和图像所有四个角的相应纬度/经度坐标.

从我们得到的Landsat元文件(ID：LC82171052016079LGN00)：

CORNER_UL_LAT_PRODUCT = -66.61490 CORNER_UL_LON_PRODUCT = -61.31816
CORNER_UR_LAT_PRODUCT = -68.74325 CORNER_UR_LON_PRODUCT = -58.04533
CORNER_LL_LAT_PRODUCT = -67.68721 CORNER_LL_LON_PRODUCT = -67.01109
CORNER_LR_LAT_PRODUCT = -69.94052 CORNER_LR_LON_PRODUCT = -64.18581
CORNER_UL_PROJECTION_X_PRODUCT = -2259300.000
CORNER_UL_PROJECTION_Y_PRODUCT = 1236000.000
CORNER_UR_PROJECTION_X_PRODUCT = -1981500.000
CORNER_UR_PROJECTION_Y_PRODUCT = 1236000.000
CORNER_LL_PROJECTION_X_PRODUCT = -2259300.000
CORNER_LL_PROJECTION_Y_PRODUCT = 958500.000
CORNER_LR_PROJECTION_X_PRODUCT = -1981500.000
CORNER_LR_PROJECTION_Y_PRODUCT = 958500.000

…

GROUP = PROJECTION_PARAMETERS MAP_PROJECTION = “PS” DATUM = “WGS84”
ELLIPSOID = “WGS84” VERTICAL_LON_FROM_POLE = 0.00000 TRUE_SCALE_LAT =
-71.00000 FALSE_EASTING = 0 FALSE_NORTHING = 0 GRID_CELL_SIZE_PANCHROMATIC = 15.00 GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00 ORIENTATION = “NORTH_UP”
RESAMPLING_OPTION = “CUBIC_CONVOLUTION” END_GROUP =
PROJECTION_PARAMETERS

通过使用具有正确地图投影的底图,我们应该能够从X / Y值导出角点lat / lon值：

import numpy as np
from mpl_toolkits.basemap import Basemap

m=Basemap(resolution='h',projection='spstere', ellps='WGS84', boundinglat=-60,lon_0=180, lat_ts=-71)

x_crn=np.array([-2259300,-1981500,-2259300,-1981500])# upper left, upper right, lower left, lower right
y_crn=np.array([1236000, 1236000, 958500, 958500])# upper left, upper right, lower left, lower right

x0, y0= m(0, -90)
#Basemap coordinates at the south pole
#note that (0,0) of the Basemap is in a corner of the map, 
#while other data sets use the south pole. 
#This is easy to take into account:
lon_crn, lat_crn = m(x0-x_crn, y0-y_crn, inverse=True)


print 'lon_crn: '+str(lon_crn)
print 'lat_crn: '+str(lat_crn)

哪个回报：

lon_crn: [-61.31816102 -58.04532791 -67.01108782 -64.1858106 ]
lat_crn: [-67.23548626 -69.3099076  -68.28071626 -70.47651326]

正如你所看到的那样,经度与图元文件中的经度一致,但纬度是低的.

我可以通过以下方式估算纬度：

lat_crn=(lat_crn+90.)*1.0275-90.

但这真的不令人满意.

如果使用图元文件中的X / Y角坐标(红色底图drawcoastlines()),则这是图像的定位方式：
《python – 为什么Basemap南极立体地图投影坐标与同一投影中的数据集不一致？》
这就是使用角落lat / lon的样子：
《python – 为什么Basemap南极立体地图投影坐标与同一投影中的数据集不一致？》

在这种情况下,我可以简单地使用纬度/经度坐标,但如前所述,有数据集(如this)仅由X / Y坐标提供,这使得依赖底图投影非常重要.我知道还有其他模块可以将数据重新投影为潜在的解决方法,但它应该可以在没有其他模块的情况下工作,并且重新投影本身可能会引入错误.

由于这个问题出现在不同的数据集中,我想相信它是底图模块中的一个错误,但我也可能会一次又一次地犯同样的错误或者有错误的期望.

最佳答案 我做了一些实验,似乎改变lat_ts对projection =’spstere’没有影响.事实上,似乎隐含地假设投影纬度为lat_ts = -90.无论你指定什么价值.

我使用projection =’stere’取得了更大的成功,因此您将在示例中构建Basemap,如下所示：

m=Basemap(width=5400000., height=5400000., projection='stere',
          ellps='WGS84', lon_0=180., lat_0=-90., lat_ts=-71.)

您可能更喜欢设置角的纬度和经度,而不是应用程序的绘图的宽度和高度.

点赞