Cartopy：海岸线不与 imshow 投影对齐
Cartopy: coastlines not lining up with imshow projection
我正在尝试在 Python 中生成一个图形，它将使用固定为正弦网格的卫星数据 (POLDER) 将 Cartopy 中的地图海岸线与 RGB 投影对齐。
Matplotlib的Basemap和Cartopy我都试过了，用Cartopy运气更好，但是即使按照别人的代码，海岸线也不匹配。
我目前拥有的：
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
proj = ccrs.Sinusoidal(central_longitude=0)
fig = plt.figure(figsize=(12, 12))
# set image extents using the min and max of lon at lat
# image_extent ~= [-73.25, -10, 59.5, 84]
min_lon = np.nanmin(subset_lon)
max_lon = np.nanmax(subset_lon)
min_lat = np.nanmin(subset_lat)
max_lat = np.nanmax(subset_lat)
extents = proj.transform_points(ccrs.Geodetic(), np.array([min_lon, max_lon]),
np.array([min_lat, max_lat]))
img_extents = [extents[0][0], extents[1][0], extents[0][1], extents[1][1]]
ax = plt.axes(projection=proj)
# # image RGB
ax.imshow(RGB, origin='upper', extent=img_extents, transform=proj))
ax.set_xmargin(0.05)
ax.set_ymargin(0.10)
ax.coastlines(color='white')
plt.show()
产生：
this figure where the coastlines do NOT match up
Figure with central_lon as -20
Figure using only imshow
我知道投影必须是正弦曲线，所以这不应该是问题。
关于它可能是什么的任何想法或关于如何修复它的提示？
这是数据集：
https://drive.google.com/file/d/1vRLKmeAXzCk5cLCJ1syOt7EJnaTmIwOa/view
以及用于提取数据并制作我想要与 cartopy 海岸线叠加的图像的代码：
data = SD(path_to_file, SDC.READ)
subset_lat = data.select('subset_lat')[:]
subset_lon = data.select('subset_lon')[:]
R = data.select('mband07')[:]
G = data.select('mband02')[:]
B = data.select('mband01')[:]
RGB = np.dstack((R, G, B))
plt.imshow(RGB)
** 编辑以添加两条注释和代码以制作 imshow RGB 图像。
谢谢！！
问题由两部分组成：
- 您需要计算投影 space 中的范围，而不是 latitude/longitude space 中的范围，因为限制坐标之间没有直接对应关系。
- 您需要知道正弦投影的参考经度，但不幸的是，它在您的数据集元数据中不可用。人工检查显示它大约在 -14 或 -15 度左右。
要通过 basemap 解决此问题，您需要一些技巧，因为目前 Basemap 不支持正弦投影的范围（即它始终设置为整个地球） .但是您可以在 Basemap 对象和基础 Axes 对象中手动更新范围，然后它应该可以工作：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from pyhdf.SD import SD
from pyhdf.SD import SDC
# Open dataset file.
path = "Greenland_PM01_L2.2007061414150564.058140.1.46-public-release.hdf"
data = SD(path, SDC.READ)
# Read latitude/longitude coordinates.
subset_lat = data.select("subset_lat")[:]
subset_lon = data.select("subset_lon")[:]
# Read RGB image.
R = data.select("mband07")[:]
G = data.select("mband02")[:]
B = data.select("mband01")[:]
RGB = np.dstack((R, G, B))
# Create `Basemap` object to plot the RGB image.
fig1 = plt.figure()
ax1 = fig1.add_axes([0, 0, 1, 1])
bmap1 = Basemap(projection="sinu", resolution="i", lon_0=-14, ax=ax1)
# Set image extents in the projected space. The hack here is that we
# update the corners by hand in both the `Axes` and `Basemap` objects.
x, y = bmap1(subset_lon, subset_lat)
llcrnrx, urcrnrx = np.nanmin(x), np.nanmax(x)
llcrnry, urcrnry = np.nanmin(y), np.nanmax(y)
bmap1.llcrnrx, bmap1.llcrnry = llcrnrx, llcrnry
bmap1.urcrnrx, bmap1.urcrnry = urcrnrx, urcrnry
ax1.set_xlim(llcrnrx, urcrnrx)
ax1.set_ylim(llcrnry, urcrnry)
# Now we can call `imshow`. The keyword argument `origin` does not work
# here, so we can simply reverse the image rows.
bmap1.drawcoastlines(color="white", linewidth=0.5)
bmap1.imshow(RGB[::-1])
# We can also play with the land cover image and check how well the
# coastlines are fitting.
landcover = data.select("subset_land")[:]
fig2 = plt.figure()
ax2 = fig2.add_axes([0, 0, 1, 1])
bmap2 = Basemap(projection="sinu", resolution="i", lon_0=-14, ax=ax2)
bmap2.llcrnrx, bmap2.llcrnry = llcrnrx, llcrnry
bmap2.urcrnrx, bmap2.urcrnry = urcrnrx, urcrnry
ax2.set_xlim(llcrnrx, urcrnrx)
ax2.set_ylim(llcrnry, urcrnry)
bmap2.drawcoastlines(color="red", linewidth=0.5)
bmap2.imshow(landcover[::-1])
显示此 RGB：
这片土地覆盖：
南方的海岸线比北方的更匹配，但我不确定这是否只是仪器本身的限制（加上正确参考经度的不确定性）。
所以我查看了数据，我认为问题实际上是因为您的数据不是在等距栅格中提供的！
因此使用 imshow 是有问题的，因为它会引入失真...（请注意，imshow 只会根据数据点的数量将您的范围划分为大小相等的像素，而不管实际坐标！）
为了阐明我的意思...我很确定您所拥有的数据在此处描述的网格中提供（第 13 页）：
https://www.theia-land.fr/wp-content/uploads/2018/12/parasol_te_level-3_format-i2.00.pdf
我个人从未遇到过这种grid-definition，但似乎调整了每个纬度的像素数以解决失真问题...引用上面的文档：
Along a parallel, the step is chosen in order to have a resolution as constant as possible. The number of pixels from 180°W to 180°E is chosen equal to 2 x NINT[3240 cos(latitude)] where NINT stands for nearest
integer.
为了快速可视化数据，我建议使用散点图 + 实际坐标而不是 imshow！
或者，我正在开发一个基于 matplotlib+ cartopy 的绘图库 (EOmaps)，旨在简化不规则或 non-uniformly 网格数据集的可视化，例如和你一样...
再投影到正弦投影后，可以看到 pixel-centers 之间的距离至少接近相等...
所以我们可以使用此信息将数据可视化为正弦投影中的矩形，大小为 6200x6200...产生预期的图像：
... 并且放大肯定会显示这些像素的大小大致相同，但它们绝对不是均匀分布的！
...这是使用 EOmaps 创建上述图的代码：
from eomaps import Maps
import numpy as np
from pyhdf.SD import SD, SDC
data = SD("Copy of greenland_PM01_L2.2007061414150564.058140.1.46-public-release.hdf", SDC.READ)
# ---- get the actual coordinates of your datapoints
lon, lat = data.select("subset_lon")[:], data.select("subset_lat")[:]
# mask nan-values
mask = np.logical_and(np.isfinite(lon), np.isfinite(lat))
lon, lat = lon[mask], lat[mask]
# ---- get the colors you want to use
R = data.select("mband07")[:][mask]
G = data.select("mband02")[:][mask]
B = data.select("mband01")[:][mask]
RGB = list(zip(R, G, B)) # a list of rgb-tuples
# ---- plot the data
m = Maps(Maps.CRS.Sinusoidal())
m.add_feature.preset.coastline()
m.set_shape.rectangles(radius=3100, radius_crs=Maps.CRS.Sinusoidal())
# use "dummy" values since we provide explicit colors
m.set_data(data=np.empty(lon.shape), xcoord=lon, ycoord=lat, crs=4326)
m.plot_map(set_extent=True, color=RGB)
...通过这种方法，海岸线也位于它们应该位于的位置！
我正在尝试在 Python 中生成一个图形，它将使用固定为正弦网格的卫星数据 (POLDER) 将 Cartopy 中的地图海岸线与 RGB 投影对齐。
Matplotlib的Basemap和Cartopy我都试过了，用Cartopy运气更好，但是即使按照别人的代码，海岸线也不匹配。
我目前拥有的：
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
proj = ccrs.Sinusoidal(central_longitude=0)
fig = plt.figure(figsize=(12, 12))
# set image extents using the min and max of lon at lat
# image_extent ~= [-73.25, -10, 59.5, 84]
min_lon = np.nanmin(subset_lon)
max_lon = np.nanmax(subset_lon)
min_lat = np.nanmin(subset_lat)
max_lat = np.nanmax(subset_lat)
extents = proj.transform_points(ccrs.Geodetic(), np.array([min_lon, max_lon]),
np.array([min_lat, max_lat]))
img_extents = [extents[0][0], extents[1][0], extents[0][1], extents[1][1]]
ax = plt.axes(projection=proj)
# # image RGB
ax.imshow(RGB, origin='upper', extent=img_extents, transform=proj))
ax.set_xmargin(0.05)
ax.set_ymargin(0.10)
ax.coastlines(color='white')
plt.show()
产生： this figure where the coastlines do NOT match up
Figure with central_lon as -20
Figure using only imshow
我知道投影必须是正弦曲线，所以这不应该是问题。
关于它可能是什么的任何想法或关于如何修复它的提示？
这是数据集： https://drive.google.com/file/d/1vRLKmeAXzCk5cLCJ1syOt7EJnaTmIwOa/view
以及用于提取数据并制作我想要与 cartopy 海岸线叠加的图像的代码：
data = SD(path_to_file, SDC.READ)
subset_lat = data.select('subset_lat')[:]
subset_lon = data.select('subset_lon')[:]
R = data.select('mband07')[:]
G = data.select('mband02')[:]
B = data.select('mband01')[:]
RGB = np.dstack((R, G, B))
plt.imshow(RGB)
** 编辑以添加两条注释和代码以制作 imshow RGB 图像。
谢谢！！
问题由两部分组成：
- 您需要计算投影 space 中的范围，而不是 latitude/longitude space 中的范围，因为限制坐标之间没有直接对应关系。
- 您需要知道正弦投影的参考经度，但不幸的是，它在您的数据集元数据中不可用。人工检查显示它大约在 -14 或 -15 度左右。
要通过 basemap 解决此问题，您需要一些技巧，因为目前 Basemap 不支持正弦投影的范围（即它始终设置为整个地球） .但是您可以在 Basemap 对象和基础 Axes 对象中手动更新范围，然后它应该可以工作：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from pyhdf.SD import SD
from pyhdf.SD import SDC
# Open dataset file.
path = "Greenland_PM01_L2.2007061414150564.058140.1.46-public-release.hdf"
data = SD(path, SDC.READ)
# Read latitude/longitude coordinates.
subset_lat = data.select("subset_lat")[:]
subset_lon = data.select("subset_lon")[:]
# Read RGB image.
R = data.select("mband07")[:]
G = data.select("mband02")[:]
B = data.select("mband01")[:]
RGB = np.dstack((R, G, B))
# Create `Basemap` object to plot the RGB image.
fig1 = plt.figure()
ax1 = fig1.add_axes([0, 0, 1, 1])
bmap1 = Basemap(projection="sinu", resolution="i", lon_0=-14, ax=ax1)
# Set image extents in the projected space. The hack here is that we
# update the corners by hand in both the `Axes` and `Basemap` objects.
x, y = bmap1(subset_lon, subset_lat)
llcrnrx, urcrnrx = np.nanmin(x), np.nanmax(x)
llcrnry, urcrnry = np.nanmin(y), np.nanmax(y)
bmap1.llcrnrx, bmap1.llcrnry = llcrnrx, llcrnry
bmap1.urcrnrx, bmap1.urcrnry = urcrnrx, urcrnry
ax1.set_xlim(llcrnrx, urcrnrx)
ax1.set_ylim(llcrnry, urcrnry)
# Now we can call `imshow`. The keyword argument `origin` does not work
# here, so we can simply reverse the image rows.
bmap1.drawcoastlines(color="white", linewidth=0.5)
bmap1.imshow(RGB[::-1])
# We can also play with the land cover image and check how well the
# coastlines are fitting.
landcover = data.select("subset_land")[:]
fig2 = plt.figure()
ax2 = fig2.add_axes([0, 0, 1, 1])
bmap2 = Basemap(projection="sinu", resolution="i", lon_0=-14, ax=ax2)
bmap2.llcrnrx, bmap2.llcrnry = llcrnrx, llcrnry
bmap2.urcrnrx, bmap2.urcrnry = urcrnrx, urcrnry
ax2.set_xlim(llcrnrx, urcrnrx)
ax2.set_ylim(llcrnry, urcrnry)
bmap2.drawcoastlines(color="red", linewidth=0.5)
bmap2.imshow(landcover[::-1])
显示此 RGB：
这片土地覆盖：
南方的海岸线比北方的更匹配，但我不确定这是否只是仪器本身的限制（加上正确参考经度的不确定性）。
所以我查看了数据，我认为问题实际上是因为您的数据不是在等距栅格中提供的！ 因此使用 imshow 是有问题的，因为它会引入失真...（请注意，imshow 只会根据数据点的数量将您的范围划分为大小相等的像素，而不管实际坐标！）
为了阐明我的意思...我很确定您所拥有的数据在此处描述的网格中提供（第 13 页）：
https://www.theia-land.fr/wp-content/uploads/2018/12/parasol_te_level-3_format-i2.00.pdf
我个人从未遇到过这种grid-definition，但似乎调整了每个纬度的像素数以解决失真问题...引用上面的文档：
Along a parallel, the step is chosen in order to have a resolution as constant as possible. The number of pixels from 180°W to 180°E is chosen equal to 2 x NINT[3240 cos(latitude)] where NINT stands for nearest integer.
为了快速可视化数据，我建议使用散点图 + 实际坐标而不是 imshow！
或者，我正在开发一个基于 matplotlib+ cartopy 的绘图库 (EOmaps)，旨在简化不规则或 non-uniformly 网格数据集的可视化，例如和你一样...
再投影到正弦投影后，可以看到 pixel-centers 之间的距离至少接近相等...
所以我们可以使用此信息将数据可视化为正弦投影中的矩形，大小为 6200x6200...产生预期的图像：
... 并且放大肯定会显示这些像素的大小大致相同，但它们绝对不是均匀分布的！
...这是使用 EOmaps 创建上述图的代码：
from eomaps import Maps
import numpy as np
from pyhdf.SD import SD, SDC
data = SD("Copy of greenland_PM01_L2.2007061414150564.058140.1.46-public-release.hdf", SDC.READ)
# ---- get the actual coordinates of your datapoints
lon, lat = data.select("subset_lon")[:], data.select("subset_lat")[:]
# mask nan-values
mask = np.logical_and(np.isfinite(lon), np.isfinite(lat))
lon, lat = lon[mask], lat[mask]
# ---- get the colors you want to use
R = data.select("mband07")[:][mask]
G = data.select("mband02")[:][mask]
B = data.select("mband01")[:][mask]
RGB = list(zip(R, G, B)) # a list of rgb-tuples
# ---- plot the data
m = Maps(Maps.CRS.Sinusoidal())
m.add_feature.preset.coastline()
m.set_shape.rectangles(radius=3100, radius_crs=Maps.CRS.Sinusoidal())
# use "dummy" values since we provide explicit colors
m.set_data(data=np.empty(lon.shape), xcoord=lon, ycoord=lat, crs=4326)
m.plot_map(set_extent=True, color=RGB)
...通过这种方法，海岸线也位于它们应该位于的位置！