当前位置：文江博客话题详情

摄影：海岸线不与Imshow投影排队

发布于 2025-01-20 00:55:38 字数 2031 浏览 3 评论 0 原文

我正在尝试在Python中产生一个数字，该人物将使用卫星数据（POLDER）将Cartopy的地图海岸线排列为固定在正弦电网上的RGB投影。

我曾尝试过Matplotlib的基部和摄影作品，并且在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()

生产：这个数字，海岸线不匹配

> ：//i.sstatic.net/bhgk6.png“ rel =“ nofollow noreferrer”> central_lon as -20

仅使用iMshow

我知道投影必须是正弦的，因此这不应该是问题。

关于它可能是什么的任何想法，或有关如何修复它的提示？

这是数据集：

和代码以提取数据并制作我想要与Cartopy Coastlines叠加的图像：

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图像。

谢谢！！

原文

I'm trying to produce a figure in Python that will line up the map coastlines from Cartopy with a RGB projection using satellite data (POLDER) that is fixed to Sinusoidal grid.

I have tried both Matplotlib's Basemap and Cartopy, and have had better luck with Cartopy, however even after following other people's code, the coastlines do not match up.

What I have so far:

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()

Produces:
this figure where the coastlines do NOT match up

Figure with central_lon as -20

Figure using only imshow

I know the projection has to be sinusoidal, so that shouldn't be the issue.

Any ideas on what else it could be or tips on how to fix it?

Here is the dataset:
https://drive.google.com/file/d/1vRLKmeAXzCk5cLCJ1syOt7EJnaTmIwOa/view

And code to extract the data and make the image that I would like overlayed with the cartopy coastlines:

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)

** Edited to add on two comments and code to make imshow RGB image.

Thanks!!

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

月寒剑心 2025-01-27 00:55:38

所以我查看了数据，我认为问题实际上来自于您的数据不是在等距栅格中提供的事实！
因此，使用 imshow 是有问题的，因为它会引入扭曲...（请注意，imshow 只会根据数据点的数量将您的范围划分为相同大小的像素，而不考虑实际坐标！）

为了澄清我的意思...我'我非常确定您所拥有的数据已在此处描述的网格中提供（第 13 页）：

https://www.theia- land.fr/wp-content/uploads/2018/12/parasol_te_level-3_format-i2.00.pdf

我个人从未遇到过这种情况网格定义，但似乎调整了每个纬度的像素数以考虑扭曲......引用上面的文档：

沿着平行线，选择步长是为了使分辨率尽可能恒定。从 180°W 到 180°E 的像素数选择等于 2 x NINT[3240 cos(latitude)]，其中 NINT 代表最近的
整数。

因此，为了快速可视化数据，我建议使用散点图+实际坐标而不是 imshow！

或者，我正在开发一个基于 matplotlib+ cartopy 的绘图库（EOmaps），旨在简化不规则或不均匀网格数据集（例如您的数据集）的可视化...

在重新投影到正弦投影后，人们可以看到像素中心之间的距离至少几乎相等...

因此，我们可以使用此信息将数据可视化为正弦投影中大小为 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)

...通过这种方法，海岸线也位于它们应该在的位置！

so I had a look at the data and I think the problem actually comes from the fact that your data is NOT provided in an equidistant raster!
Therefore using imshow is problematic since it will introduce distortions... (note that imshow will just divide your extent into equally sized pixels based on the number of datapoints, irrespective of the actual coordinates!)

To clarify what I mean... I'm pretty sure the data you have is provided in the grid described here (page 13):

https://www.theia-land.fr/wp-content/uploads/2018/12/parasol_te_level-3_format-i2.00.pdf

I have personally never encountered this kind of grid-definition, but it seems that the number of pixels per latitude is adjusted to account for distortions... to quote from the document above:

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.

To quickly visualize the data I'd therefore suggest to use a scatterplot + the actual coordinates instead of imshow!

Alternatively, I'm developing a matplotlib+ cartopy based plotting library (EOmaps) that is intended to simplify visualizations from irregular or non-uniformly gridded datasets such as yours...

After reprojection to Sinusoidal projection, one can see that the distance between the pixel-centers is at least nearly equal...

So we can use this information to visualize the data as rectangles in Sinusoidal projection with a size of 6200x6200... yielding the expected image:

... and zooming in certainly shows that those pixels have approximately the same size but they are definitely not equally distributed!

... here's the code to create the above plot with 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)

... and with this approach also the coastlines are located where they are supposed to be!

回复收藏 0 原文

债姬 2025-01-27 00:55:38

该问题由两部分组成：

您需要计算投影空间中的范围，而不是纬度/经度空间中的范围，因为限制坐标之间没有直接对应关系。
您需要知道正弦投影的参考经度，但不幸的是它在数据集元数据中不可用。手动检查显示温度约为 -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：

和土地覆盖：

南部的海岸线比北部的匹配得更好，但我不确定这是否只是仪器本身的限制（加上正确参考经度的不确定性）。

The issue consists of two parts:

You need to compute the extents in the projected space and not in the latitude/longitude space, because there is not a direct correspondence between limit coordinates.
You need to know the reference longitude of your sinusoidal projection, but unfortunately it is not available in your dataset metadata. A manual inspection showed me that it is somewhere around -14 or -15 degrees.

To solve this by means of basemap, you need a bit of hacking because at the moment Basemap does not support extents for sinusoidal projection (i.e. it is always set to the whole globe). But you can update the extent by hand both in the Basemap object and the underlying Axes object and then it should work:

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])

which shows this RGB:

and this land cover:

The coastlines match better in the South than in the North, but I am not sure if this is just a limitation of the instrument itself (plus the uncertainty in the correct reference longitude).

回复收藏 0 原文

~没有更多了~