Блог

Главная - Вопросы по программированию - Неправильная пространственная информация после повторной выборки 2D пространственных данных (разного происхождения и разрешения) с помощью pyresample в Python

Неправильная пространственная информация после повторной выборки 2D пространственных данных (разного происхождения и разрешения) с помощью pyresample в Python

#python-3.x #matplotlib #geospatial #resampling #cartopy

#python-3.x #matplotlib #геопространственный #повторная выборка #картография

Вопрос:

Мне нужно регулярно выполнять повторную выборку данных с сеткой (в lon-lat) в новую сетку с более низким разрешением и другим источником. Я бы использовал 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, 0.4, 0.6, 0.5, 1., 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, xmax, ymin, 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,2,1,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', proj_id='proj_id_blabla',
                            projection='EPSG:4326', 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, 1, 1, 0, 0)], [targetextentforplot[y] for y in (2, 2, 3, 3, 2)],
            'r--', lw=3, transform=ccrs.PlateCarree())
    ax.set_extent([-12, 12, 38, 62])

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

doplt(ax1, sourcedata, 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 , в идеале с example.

Ответ №1:

Итак, проблема в том, что вы предполагаете, что строка 0 является нижней частью изображения, но, как показано в этом примере, 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, 0.4, 0.6, 0.5, 1., 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, xmax, ymin, 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,2,1,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', proj_id='proj_id_blabla',
                            projection='EPSG:4326', 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, 1, 1, 0, 0)], [targetextentforplot[y] for y in (2, 2, 3, 3, 2)],
            'r--', lw=3, transform=ccrs.PlateCarree())
    ax.set_extent([-12, 12, 38, 62])

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

doplt(ax1, sourcedata, 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?)');

Это дает:

Результирующее изображение

Я счел полезным использовать изображение с большим разнообразием, чтобы иметь возможность лучше сопоставлять его с исходными данными.

Комментарии:

1. Спасибо! Я проверю это. Но, делая вид, что данные примера являются пространственными данными, результат здесь будет неправильным (пространственная информация остается потерянной). Итак, предполагается ли, что в таком случае мы должны фактически переместить их по широте, прежде чем что-либо делать в pyresample , а затем преобразовать их снова?! Возможно, вы можете изменить свой ответ, чтобы он точно соответствовал пространственной информации из вопроса (как если бы это была карта, здесь страны были бы неправильными)?

2. Я не уверен, что вы подразумеваете под «пространственная информация остается потерянной». Единственное изменение, которое я внес, заключалось в том, что строка 0 — это высокая широта. Это расположение точно так, как это сделано в исходном примере: pyresample.readthedocs.io/en/latest/… Это также то, что указано для ImageDraw : pillow.readthedocs.io/en/stable/reference/… Действительно, проблема здесь заключалась в том, как вы вычисляли широту по координатам полигона (которые были в долях изображения) Часть 0 должна отображаться на latmax , а не latmin .

3. Спасибо — извините, я не был яснее. Я имел в виду, что если мой исходный массив был фактическими пространственными данными, ваша манипуляция подразумевает «вертикальный флип» (т. Е. Продолжать охватывать в основном Францию, как изначально, а не Северное море, как в вашем решении), но это просто операция флип, по крайней мере, resample операция работает, и вы заметили проблему — спасибо! [Тем временем я обнаружил, что использование SwathDefinition немедленно решило проблему …]

Ответ №2:

Я обнаружил, что использование SwathDefinition вместо AreaDefinition (см. Документ) решило проблему.

Определение 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)

введите описание изображения здесь