cartopy：NOAA APT 图像上的地图叠加层

Question

我正在进行一个尝试解码 NOAA APT 图像的项目，到目前为止，我已经到了可以从 RTLSDR 的原始 IQ 记录中获取图像的阶段。这是解码后的图像之一， Decoded NOAA APT image此图像将用作代码的输入（此处显示为 m3.png）

现在我正在努力在图像上覆盖地图边界（注意：仅在上图的左半部分）

我们知道，图像被捕获的时间和卫星信息：位置、方向等。所以，我使用卫星的位置来获取地图投影的中心和卫星的方向来适当地旋转图像.

首先我在 Basemap 中尝试过，这里是代码

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
from scipy import ndimage

im = plt.imread('m3.png')
im = im[:,85:995] # crop only the first part of whole image

rot = 198.3913296679117 # degrees, direction of sat movement
center = (50.83550180700588, 16.430852851867176) # lat long

rotated_img = ndimage.rotate(im, rot) # rotate image

w = rotated_img.shape[1]*4000*0.81 # in meters, spec says 4km per pixel, but I had to make it 81% less to get better image
h = rotated_img.shape[0]*4000*0.81 # in meters, spec says 4km per pixel, but I had to make it 81% less to get better image

m = Basemap(projection='cass',lon_0 = center[1],lat_0 = center[0],width = w,height = h, resolution = "i")
m.drawcoastlines(color='yellow')
m.drawcountries(color='yellow')

im = plt.imshow(rotated_img, cmap='gray', extent=(*plt.xlim(), *plt.ylim()))
plt.show()

结果我得到了this image，看起来还不错

我想将代码移至 Cartopy，因为它更易于安装并且正在积极开发中。我找不到类似的方法来设置边界，即以米为单位的宽度和高度。所以，我修改了most similar example。我找到了一个函数，可以将米添加到经纬度，并用它来设置边界。

这是 Cartopy 中的代码，

import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from scipy import ndimage
import cartopy.feature

im = plt.imread('m3.png')
im = im[:,85:995] # crop only the first part of whole image

rot = 198.3913296679117 # degrees, direction of sat movement
center = (50.83550180700588, 16.430852851867176) # lat long

def add_m(center, dx, dy):
    # source: https://stackoverflow.com/questions/7477003/calculating-new-longitude-latitude-from-old-n-meters
    new_latitude  = center[0] + (dy / 6371000.0) * (180 / np.pi)
    new_longitude = center[1] + (dx / 6371000.0) * (180 / np.pi) / np.cos(center[0] * np.pi/180)
    return [new_latitude, new_longitude]

fig = plt.figure()

img = ndimage.rotate(im, rot)

dx = img.shape[0]*4000/2*0.81 # in meters
dy = img.shape[1]*4000/2*0.81 # in meters

leftbot = add_m(center, -1*dx, -1*dy)
righttop = add_m(center, dx, dy)

img_extent = (leftbot[1], righttop[1], leftbot[0], righttop[0])

ax = plt.axes(projection=ccrs.PlateCarree())
ax.imshow(img, origin='upper', cmap='gray', extent=img_extent, transform=ccrs.PlateCarree())
ax.coastlines(resolution='50m', color='yellow', linewidth=1)
ax.add_feature(cartopy.feature.BORDERS, linestyle='-', edgecolor='yellow')

plt.show()

Here is the result来自 Cartopy，不如 Basemap 的结果。

我有以下问题：

1. 我发现无法旋转地图而不是图像，在 底图和cartopy。因此我求助于旋转图像，是 有没有办法旋转地图？
2. 如何提高cartopy的输出？我认为这是进入的方式 我正在计算问题的程度。有没有办法我可以 提供米来设置图像的边界？
3. 有没有更好的方法来做我想做的事情？任何特定于此类应用的投影？
4. 我正在手动调整比例（我决定每像素公里数的部分），有没有办法根据 在 卫星的高度？

任何形式的输入都将受到高度赞赏。非常感谢您的宝贵时间！

如果你有兴趣可以找到项目here。

Answer 1

据我所知，底层 Proj.4 无法定义具有旋转视角的卫星投影（很高兴以其他方式显示 - 我不是专家！）（注意：也许通过 ob_tran? ）。这是您无法使用 Basemap 或 Cartopy 在“本机”坐标/方向中执行此操作的主要原因。

这个问题实际上归结为一个地理配准问题，我在 https://www.cder.dz/download/Art7-1_1.pdf 这样的地方找不到足够的信息。

我的解决方案完全是胡扯，但确实让您非常接近引用 this 图像。我加倍的软糖因素实际上是通用的，如果您想编写通用代码，这有点问题。

我不得不做的一些软糖（反复试验）：

• 将卫星方位调整 3.2 度
• 通过将图像中心沿卫星轨迹移动 10 公里来调整图像中心的位置
• 通过沿卫星轨迹垂直移动 10 公里来调整图像中心的位置
• 将 x 和 y 像素大小分别缩放 0.62 和 0.65
• 在不切实际的卫星高度上使用“近侧透视”投影

结果似乎是一个相对良好的注册图像，但正如我所说，似乎不太可能普遍适用于收到的所有图像：

生成此图像的代码（相当复杂，但完整）：

import urllib.request
urllib.request.urlretrieve('https://i.stack.imgur.com/UBIuA.jpg', 'm3.jpg')


import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from scipy import ndimage
import cartopy.feature

im = plt.imread('m3.jpg')
im = im[:,85:995] # crop only the first part of whole image

rot = 198.3913296679117 # degrees, direction of sat movement
center = (50.83550180700588, 16.430852851867176) # lat long

import numpy as np
from cartopy.geodesic import Geodesic
import matplotlib.transforms as mtransforms
from matplotlib.axes import Axes


tweaked_rot = rot - 3.2

geod = Geodesic()
# Move the center along the trajectory of the satellite by 10KM
f = np.array(
        geod.direct([center[1], center[0]],
                    180 - tweaked_rot,
                    10000))

tweaked_center = f[0, 0], f[0, 1]

# Move the satellite perpendicular from its proposed trajectory by 15KM
f = np.array(
        geod.direct([tweaked_center[0], tweaked_center[1]],
                    180 - tweaked_rot + 90,
                    10000))
tweaked_center = f[0, 0], f[0, 1]


data_crs = ccrs.NearsidePerspective(
    central_latitude=tweaked_center[1],
    central_longitude=tweaked_center[0],
)

# Compute the center in data_crs coordinates.
center_lon_lat_ortho = data_crs.transform_point(
    tweaked_center[0], tweaked_center[1], ccrs.Geodetic())

# Define the affine rotation in terms of matplotlib transforms.
rotation = mtransforms.Affine2D().rotate_deg_around(
    center_lon_lat_ortho[0], center_lon_lat_ortho[1], tweaked_rot)

# Some fudge factors. Sorry - there are entirely application specific,
# perhaps some reading of https://www.cder.dz/download/Art7-1_1.pdf
# would enlighten these... :(
ff_x, ff_y = 0.62, 0.65
ff_x = ff_y = 0.81
x_extent = im.shape[1]*4000/2 * ff_x
y_extent = im.shape[0]*4000/2 * ff_y
img_extent = [-x_extent, x_extent, -y_extent, y_extent]


fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=data_crs)
ax.margins(0.02)
with ax.hold_limits():
    ax.stock_img()


# Uing matplotlib's image transforms if the projection is the
# same as the map, otherwise we need to fall back to cartopy's
# (slower) image resampling algorithm
if ax.projection == data_crs:
    transform = rotation + ax.transData
else:
    transform = rotation + data_crs._as_mpl_transform(ax)


# Use the original Axes method rather than cartopy's GeoAxes.imshow.
mimg = Axes.imshow(ax, im, origin='upper', cmap='gray',
                   extent=img_extent, transform=transform)


lower_left = rotation.frozen().transform_point([-x_extent, -y_extent])
lower_right = rotation.frozen().transform_point([x_extent, -y_extent])
upper_left = rotation.frozen().transform_point([-x_extent, y_extent])
upper_right = rotation.frozen().transform_point([x_extent, y_extent])

plt.plot(lower_left[0], lower_left[1],
         upper_left[0], upper_left[1],
         upper_right[0], upper_right[1],
         lower_right[0], lower_right[1],
         marker='x', color='black',
         transform=data_crs)

ax.coastlines(resolution='10m', color='yellow', linewidth=1)
ax.add_feature(cartopy.feature.BORDERS, linestyle='-', edgecolor='yellow')


sat_pos = np.array(geod.direct(tweaked_center, 180 - tweaked_rot,
                               np.linspace(-x_extent*2, x_extent*2, 50)))

with ax.hold_limits():
    plt.plot(sat_pos[:, 0], sat_pos[:, 1], transform=ccrs.Geodetic(),
             label='Satellite path')
plt.plot(tweaked_center, 'ob')
plt.legend()

您可能会说，我对这个问题有点忘乎所以。这是一个非常有趣的问题，但并不是一个真正的cartopy/Basemap。

希望有帮助！