使用pyresample重新采样并绘制从geos到墨卡托的satpy场景数据

Question

我有一个来自 eumetsat 卫星的场景对象。

我感兴趣的领域： （xmin_corner、ymin_corner、xmax_corner、ymax_corner） [40,-10,105,40] 度

看起来像这样

区域定义如下所示：

Area ID: msg_seviri_unknown_3km
Description: MSG SEVIRI unknown area definition with 3 km resolution
Projection: {'a': '6378169', 'h': '35785831', 'lon_0': '45.5', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1848
Number of rows: 1665
Area extent: (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914)

我想将其投影到墨卡托投影中以进行平面观看。

我使用上面的场景区域详细信息创建一个区域定义。

                my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
                            width=1848, height=1665,
                            area_extent= (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914))
                
                new_scn = scn.resample(my_area)
                
                new_scn.show('natural_color')

并检查新数据。

投影看起来不错，但某些区域被裁剪了。

然后我尝试明确指定区域范围度数

my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
            width=1848, height=1665,
            area_extent= [40,-10,105,40], units='degrees')

new_scn = scn.resample(my_area)

new_scn.show('natural_color')

现在我得到了更好的图像，但有一些倾斜。它在水平轴上收缩。 我知道会发生这种情况，因为我已经指定了宽度和高度？

如何选择最佳高度、宽度以获得最佳真实输出投影后效果？

现在，我想生成一个图像，当我从 geos 移动到 Merc 时保持逼真的外观（长宽比？），而不指定一些高度和宽度并倾斜图像。 我希望图像看起来像这样

附加信息： 我认为因为我的数据集有 np.nans，因为它位于卫星视图的边缘，所以这会产生一些问题。

当我将区域范围从原始场景转换为经纬度时。

scn_area = scn['natural_color'].attrs['area']
# xmin, ymin
scn_area.get_lonlat_from_projection_coordinates(-604581.2379, -1096647.3571)
# (39.95708481758364, -10.007507617303142) - same as original. 40, -10

# xmax, ymax
scn_area.get_lonlat_from_projection_coordinates(-604581.2379, -1096647.3571)
# (inf, inf)

这是计算区域范围时的预期行为吗？我不应该从投影坐标中明确得到 lon lat ，而不是 inf 吗？

我猜这个 inf inf 是因为拐角处的 nan 值，在 geos 投影中它不是地球，而是空间。

现在。 当我对新的 Merc 区域定义执行相同操作时，我得到了这个。

my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
            width=1848, height=1665,
            area_extent= (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914))

#xmin, ymin
my_area.get_lonlat_from_projection_coordinates(-604581.2379, -1096647.3571)
# (40.068954335025296, -9.867939236404165) - close to original (40, -10)

#xmax, ymax
my_area.get_lonlat_from_projection_coordinates(4940163.8125, 3899023.914)
# (89.87824658822916, 33.2040663665062) - completely off. and im guessing this is why the clipping happens. 

我意识到这一定与初始化 Scene 对象时自动计算区域范围或其他原因有关，因为外部空间 np.nans。 通过明确提及我的面积范围（以度为单位），我能够克服这个问题。

现在。假设我所做的和我的理解是正确的。 我希望能够获得原始长宽比的最终图像（如果我可以这样称呼它），或者谷歌地图等上使用的最佳实践是什么。

Answer 1

这里有很多不同的变量在起作用，您已经在一个问题中涵盖了相当多的场景。我会尽力回答。

首先，我不清楚您是否试图将某些外部地理区域定义与您的区域定义相匹配。因此，如果这里有您不能/不想更改的变量，请原谅我。在第一个以米为单位定义范围的情况下，您说有些东西被“剪裁”了。为何如此？如果是的话，那为什么不增加你的范围呢？您所在区域的分辨率为 ~3km x ~3km。您可以从调用中删除

width/height

并指定

resolution=(3000, 3000)

，然后修改范围并让 Pyresample 计算出需要多少像素。

在第二种情况下，度数范围会变得棘手。您可以打印该区域和

pixel_size_x/y

属性，可以看到像素分辨率比您原来的 3km 区域大得多，并且它们在 x 和 y 维度上不相等：

In [6]: my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
   ...:             width=1848, height=1665,
   ...:             area_extent= [40,-10,105,40], units='degrees')

In [7]: my_area
Out[7]:
Area ID: my_area
Description: my_area
Projection: {'datum': 'WGS84', 'k': '1', 'lon_0': '45.5', 'no_defs': 'None', 'proj': 'merc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1848
Number of rows: 1665
Area extent: (-612257.1994, -1111475.1029, 6623509.7022, 4838471.3981)

In [8]: my_area.pixel_size_x
Out[8]: 3915.458280066441

In [9]: my_area.pixel_size_y
Out[9]: 3573.5414419900067

右侧奇怪的缺失数据弧可能是由于 Satpy 中的“reduce_data”功能造成的。这是我们正在努力改进的一个持续存在的问题。您可以通过将

reduce_data=False

传递给

Scene.resample

调用来解决这个问题，看看它是否会改善这一点。如果这不能解决问题，请告诉我。

在你的最后一组问题中，你说墨卡托投影坐标到经度/纬度“完全关闭”。你是什么意思？为何如此？你是对的，对地静止区域的 Infs 来自空间像素。一般来说，您必须注意您使用的单位以及这些坐标在什么坐标系中有效。对地静止投影中的米与墨卡托投影中的米不同。转换为经度/纬度只会进一步混淆事情，因为对地静止投影中的最小/最大 x/y 与度空间中的最小/最大经度/纬度无关，并且可能不会是墨卡托空间中的最小/最大。由于它们是不同的坐标系，因此您无法保证一个坐标系中图像的最左侧部分是另一个坐标系中图像的最左侧部分。希望这些漫无目的的内容有一定道理。

我的建议：如果您知道您想要墨卡托中的 3km x 3km 区域并希望它覆盖度空间中的特定区域，那么：

import xarray as xr
my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
    resolution=xr.DataArray([3000, 3000], attrs={"units": "meters"}),
    area_extent= [40,-10,105,40], units='degrees')

我认为这应该有效。它基本上是说“墨卡托投影中的范围以度为单位，分辨率以米为单位”。

我的猜测是

reduce_data=False

还需要解决裁剪/裁剪地球静止数据的问题。