如何解决在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)
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。