我正在进行一个尝试解码 NOAA APT 图像的项目，到目前为止，我已经到了可以从 RTLSDR 的原始 IQ 记录中获取图像的阶段。这是解码后的图像之一，已 解码的 NOAA APT 图像，此图像将用作代码的输入（此处显示为 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()
结果我得到了这张图片，看起来还不错
我想将代码移至 Cartopy，因为它更易于安装并且正在积极开发中。我找不到类似的方法来设置边界，即以米为单位的宽度和高度。所以，我修改了最相似的例子。我找到了一个函数，可以在 longs 和 lats 中添加米，并用它来设置边界。
这是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()
这是Cartopy 的结果，它不如 Basemap 的结果。
我有以下问题：
- 我发现在底图和 cartopy 中都无法旋转地图而不是图像。因此我求助于旋转图像，有没有办法旋转地图？
- 如何提高cartopy的输出？我认为这是我计算问题程度的方式。有没有办法可以提供米来设置图像的边界？
- 有没有更好的方法来做我想做的事情？任何特定于此类应用程序的投影？
- 我正在手动调整比例（我决定每像素公里数的部分），有没有办法根据卫星的高度来做到这一点？
任何形式的输入将不胜感激。非常感谢您的参与！
如果您有兴趣，可以在这里找到该项目。
