在Python中使用pyresample重新采样2D空间数据不同的原点和分辨率后，空间信息错误

问题描述

我需要将定期网格化的数据（经纬度）重新采样到分辨率较低且来源不同的新网格。我虽然会使用pyresample。

问题：重新采样后，我得到的结果的空间位置明显错误。 在以下示例中，我构造了一个简单的2D数组，其中包含一些空间网格（在sourcegrid中定义为pyresample AreaDeFinition对象）和一些遮罩，以将其重新采样到另一个targetgrid。空间信息在此过程中丢失了某个地方，我不知道在哪里...任何想法？

import numpy as np
from pyresample.geometry import AreaDeFinition
from pyresample.kd_tree import resample_nearest
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from PIL import Image,ImageDraw

# Source data
lonmin = -10; lonmax = 10.; latmin=40.; latmax=60.; nlon = 300; nlat = 250
lon = np.linspace(lonmin,lonmax,nlon); lat = np.linspace(latmin,latmax,nlat)
dlon = lon[1] - lon[0]; dlat = lat[1] - lat[0]
lon2d,lat2d = np.meshgrid(lon,lat)
sourcedata = np.cos(np.deg2rad(lat2d)*100) + np.sin(np.deg2rad(lon2d)*100)

# Introduce a polygon as mask
xpol = [frac*(nlon-1) for frac in (0,0.5,0.4,0.6,0.9,0.,0)]
ypol = [frac*(nlat-1) for frac in (0,1.,0)]
polygon = [xy for xy in zip(xpol,ypol)]
img = Image.new('L',(nlon,nlat),0)
ImageDraw.Draw(img).polygon(polygon,outline=1,fill=1)
mask = np.array(img)
xpol = [lon[int(x)] for x in xpol]; ypol = [lat[int(y)] for y in ypol] # translate in lon-lat for plot
sourcedata = np.ma.masked_where(mask,sourcedata)


# Define source and target areas
sourceextent = [lonmin-dlon/2,latmin-dlat/2,lonmax+dlon/2,latmax+dlat/2] # [xmin,ymin,xmax,ymax]
sourceextentforplot = [sourceextent[i] for i in (0,2,1,3)] # [xmin,ymax]
targetextent =  [lonmin-dlon/2 + 0.12*(lonmax-lonmin),latmin-dlat/2 + 0.24*(latmax-latmin),lonmin-dlon/2 + 0.78*(lonmax-lonmin),latmin-dlat/2 + 0.91*(latmax-latmin)]
targetextentforplot = [targetextent[i] for i in (0,3)] 

sourcegrid = AreaDeFinition(area_id='Grd1',description='Source Grid',proj_id='proj_id_blabla',projection='epsg:4326',width=nlon,height=nlat,area_extent=sourceextent)
# Lower resolution,different origin
targetgrid = AreaDeFinition(area_id='Grd2',description='Target Grid',width=123,height=97,area_extent=targetextent)

# Resample sourcedata to newdata
newdata = resample_nearest(sourcegrid,sourcedata,targetgrid,fill_value=None,radius_of_influence=50000)

# Plot
def doplt(ax,data,extent):
    ax.coastlines(resolution='50m',color='gray',alpha=1.,linewidth=2.)
    ax.gridlines(draw_labels=True)
    ax.imshow(data,origin='lower',transform=ccrs.PlateCarree(),extent=extent)
    ax.plot(xpol,ypol,'k--',transform=ccrs.PlateCarree())
    ax.plot([targetextentforplot[x] for x in (0,0)],[targetextentforplot[y] for y in (2,3,2)],'r--',lw=3,transform=ccrs.PlateCarree())
    ax.set_extent([-12,12,38,62])

fig,(ax1,ax2) = plt.subplots(2,figsize=(5,10),subplot_kw={'projection': ccrs.PlateCarree()})

doplt(ax1,extent=sourceextentforplot)
ax1.set_title('Source data,target area in red')
doplt(ax2,newdata,extent=targetextentforplot)
ax2.set_title('New data,with wrong spatial ref (or plotting?)')

plt.show()

注意：欢迎您提出除pyresample以外的其他建议来进行重采样操作。

解决方法

因此，问题在于您假设第0行位于图像的底部，但是如this example所示，pyresample使用第0行作为顶部。我修改了您的示例以调整多边形的纬度，并使用origin='upper'与imshow进行绘制：

import numpy as np
from pyresample.geometry import AreaDefinition
from pyresample.kd_tree import resample_nearest
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from PIL import Image,ImageDraw

# Source data
lonmin = -10; lonmax = 10.; latmin=40.; latmax=60.; nlon = 300; nlat = 250
lon = np.linspace(lonmin,lonmax,nlon); lat = np.linspace(latmin,latmax,nlat)
dlon = lon[1] - lon[0]; dlat = lat[1] - lat[0]
lon2d,lat2d = np.meshgrid(lon,lat)
sourcedata = np.hypot(lon2d,lat2d - 50) * (np.cos(np.deg2rad(lat2d)*100) + np.sin(np.deg2rad(lon2d)*100))

# Introduce a polygon as mask
xpol = [frac*(nlon-1) for frac in (0,0.5,0.4,0.6,0.9,0.,0)]
ypol = [frac*(nlat-1) for frac in (0,1.,0)]
polygon = [xy for xy in zip(xpol,ypol)]
img = Image.new('L',(nlon,nlat),0)
ImageDraw.Draw(img).polygon(polygon,outline=1,fill=1)
mask = np.array(img)
xpol = [lon[int(x)] for x in xpol]; ypol = [lat[nlat - 1 - int(y)] for y in ypol] # translate in lon-lat for plot
sourcedata = np.ma.masked_where(mask,sourcedata)


# Define source and target areas
sourceextent = [lonmin-dlon/2,latmin-dlat/2,lonmax+dlon/2,latmax+dlat/2] # [xmin,ymin,xmax,ymax]
sourceextentforplot = [sourceextent[i] for i in (0,2,1,3)] # [xmin,ymax]
targetextent =  [lonmin-dlon/2 + 0.12*(lonmax-lonmin),latmin-dlat/2 + 0.24*(latmax-latmin),lonmin-dlon/2 + 0.78*(lonmax-lonmin),latmin-dlat/2 + 0.91*(latmax-latmin)]
targetextentforplot = [targetextent[i] for i in (0,3)] 

sourcegrid = AreaDefinition(area_id='Grd1',description='Source Grid',proj_id='proj_id_blabla',projection='EPSG:4326',width=nlon,height=nlat,area_extent=sourceextent)
# Lower resolution,different origin
targetgrid = AreaDefinition(area_id='Grd2',description='Target Grid',width=123,height=97,area_extent=targetextent)

# Resample sourcedata to newdata
newdata = resample_nearest(sourcegrid,sourcedata,targetgrid,fill_value=None,radius_of_influence=50000)

# Plot
def doplt(ax,data,extent):
    ax.coastlines(resolution='50m',color='gray',alpha=1.,linewidth=2.)
    ax.gridlines(draw_labels=True)
    ax.imshow(data,transform=ccrs.PlateCarree(),extent=extent,norm=plt.Normalize(0,20),origin='upper')
    ax.plot(xpol,ypol,'k--',transform=ccrs.PlateCarree())
    ax.plot([targetextentforplot[x] for x in (0,0)],[targetextentforplot[y] for y in (2,3,2)],'r--',lw=3,transform=ccrs.PlateCarree())
    ax.set_extent([-12,12,38,62])

fig,(ax1,ax2) = plt.subplots(2,figsize=(5,10),subplot_kw={'projection': ccrs.PlateCarree()})

doplt(ax1,extent=sourceextentforplot)
ax1.set_title('Source data,target area in red')
doplt(ax2,newdata,extent=targetextentforplot)
ax2.set_title('New data,with wrong spatial ref (or plotting?)');

给出：

我发现使用变化较大的图像以使其与源数据更好地对齐很有帮助。

,

我发现使用SwathDefinition代替AreaDefinition（请参阅doc）解决了这个问题。

在原始代码中如下定义sourcegrid和targetgrid会得到很好的结果：

sourcegrid = SwathDefinition(lons=lon2d,lats=lat2d)
lon2dtarget,lat2dtarget = np.meshgrid(np.linspace(targetextent[0],targetextent[2],123),np.linspace(targetextent[1],targetextent[3],97))
targetgrid = SwathDefinition(lons=lon2dtarget,lats=lat2dtarget)