pytroll / pyresample

Geospatial image resampling in Python
http://pyresample.readthedocs.org
GNU Lesser General Public License v3.0
351 stars 94 forks source link

Unexpected results resampling Lambert Conformal to PlateCarree: pyresample or cartopy problem? #416

Closed hgordo closed 2 years ago

hgordo commented 2 years ago

Code Sample, a minimal, complete, and verifiable piece of code

# Your code here
import numpy as np
from pyresample import geometry,area_config
from pyresample.kd_tree import resample_nearest
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LogNorm
#from netCDF4 import Dataset
proj_string='+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
lower_left_x = -2412000
lower_left_y = -1620000
upper_right_x = -2412000+12000*396
upper_right_y = -1620000+12000*246
msg_area = geometry.AreaDefinition('epa_12us2', 'EPA 12 US 2 grid','epa_12us2',proj_string,396, 246,[lower_left_x,lower_left_y,upper_right_x,upper_right_y])

reff =np.ones((396,246))
#d = Dataset('emis_mole_all_20140802_12US2_nobeis_2014fd_nata_cb6_14j_nohap.nc4')
#reff= np.array(d['NH3'])[12,0,:,:]
area_def = geometry.AreaDefinition.from_extent(area_id='conus',projection='+proj=eqc +ellps=WGS84',shape=(400,300),area_extent=(-140,20,-60,60),units='deg')
result_data = resample_nearest(msg_area,reff,area_def,fill_value=-99,radius_of_influence=50000)

crs = msg_area.to_cartopy_crs()
newcrs=area_def.to_cartopy_crs()
print(crs)

fig, ax = plt.subplots(subplot_kw=dict(projection=crs))
plt.imshow(reff, transform = crs,extent=crs.bounds,origin='lower',norm=LogNorm())
plt.colorbar()
ax.coastlines()
ax.gridlines()

fig, ax = plt.subplots(subplot_kw=dict(projection=newcrs))
plt.imshow(result_data, transform = newcrs,extent=newcrs.bounds,origin='lower',norm=LogNorm())
plt.gca().coastlines()
ax.gridlines()
cbar = plt.colorbar()
plt.show()

Problem description

When resampling the area definitions in the above, I see image results from image

This is not obviously wrong, but I am axctually trying to regrid an inventory of ammonia emissions, which look like image before regridding and image after. I have included the code to read the data file, commented out, in my sample above. The file is available online but it is 860MB.

I am not sure if I am mis-using the software - apologies!

Expected Output

I expect the green area of 1 to cover Florida in the resampled output

Actual Result, Traceback if applicable

The resampled green area seems distorted, it doesn't cover Florida for example, if the Cartopy coastlines are right.

Versions of Python, package at hand and relevant dependencies

python 3.9.9 h62f1059_0_cpython conda-forge pyresample 1.22.3 py39hde0f152_0 conda-forge

hgordo commented 2 years ago

i think I am misinterpreting what pyresample can do, sorry for the noise on your github page! Thanks.

djhoese commented 2 years ago

Based on your code I don't think there is anything you are doing that is out of the question. I'd be OK with reopening this (or just continuing the discussion with it closed).

I don't see anything super wrong in the code. It all looks correct at first glance. However, I'm wondering if your area description of your input is correct. You have the images from when you resampled the area definitions, how did you generate that? You can see the resampled area doesn't even cover Florida. But the not resampled data looks fine... :thinking: this is very strange. If you can provide more of the code for how you plotted the areas it would be appreciated.

djhoese commented 2 years ago

Ah sorry, I just realized your code is the plotting of the grids with np.ones.

djhoese commented 2 years ago

I got it! Or at least I got something that made sense to me. Change your origin to origin="upper". That fixed your code for me above. If your image data ends up looking wrong with this, then you may need to flip it in the Y direction data = data[::-1, :] before resampling it. You could also flip the source area definition's Y extents, that should work too. My guess is your image data is flipped from what pyresample assumes which is typically that the first pixel of an array is the upper-left corner of the image.

hgordo commented 2 years ago

hi, Indeed that fixed it! Thank you very much for the help and sorry again for bothering you. I hadn't figured out that the image data were indeed flipped in the Y direction. Clearly pyresample is now doing a great job.

Screenshot 2022-02-02 at 09 00 59