Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
hgordo opened this issue Feb 1, 2022 · 5 comments

Comments

@hgordo
Copy link

hgordo commented Feb 1, 2022

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
Copy link
Author

hgordo commented Feb 2, 2022

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

@hgordo hgordo closed this as completed Feb 2, 2022
@djhoese
Copy link
Member

djhoese commented Feb 2, 2022

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... 🤔 this is very strange. If you can provide more of the code for how you plotted the areas it would be appreciated.

@djhoese
Copy link
Member

djhoese commented Feb 2, 2022

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

@djhoese
Copy link
Member

djhoese commented Feb 2, 2022

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
Copy link
Author

hgordo commented Feb 2, 2022

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants