resample from latlon large target grid to smaller grid

[Marston]

I'm trying to resample and reproject some data from a global, regular latlon grid to much smaller LAEA grid, which is essentially going from a global gris defined in degrees to a smaller region defined in meters. The problem comes in here:

result = kd_tree.resample_nearest(grid_def, da.values, area_def, radius_of_influence=50000, epsilon=0.5)

Code example

import os
os.environ["PROJ_LIB"] = "/root/anaconda3/envs/exp/share/proj"
import json
import xarray as xr
import numpy as np
from pyresample.geometry import AreaDefinition, GridDefinition
from pyresample import kd_tree, save_quicklook

def prj4_string_to_dict(pj4str):
    ret = {}
    segs = pj4str.split(' ')
    for s in segs:
        if s == '': continue
        if '+units' in s:
            ret['units'] = s.split('=')[-1]
        elif '+proj' in s:
            ret['proj'] = s.split('=')[-1]
        elif '+lat_0' in s:
            ret['lat_0'] = s.split('=')[-1]
        elif '+lon_0' in s:
            ret['lon_0'] = s.split('=')[-1]
        elif '+a' in s:
            ret['R'] = s.split('=')[-1]
        elif '+ellps' in s:
            ret['ellps'] = s.split('=')[-1]
        elif '+over' in s:
            ret['over'] = s.split('=')[-1]
        else:
            print(f'No logic for {s}')
            print(segs)
    return ret
# Load the data array from the netCDF file
da = xr.open_dataarray('gfs_orography.nc', engine='netcdf4')
lon2d, lat2d = np.meshgrid(da.coords['lon'].values, da.coords['lat'].values)
# there it is required that the coordinates arrays are the same shape, which necessitates the following: 
grid_def = GridDefinition(lons=lon2d, lats=lat2d)

# Load the region projection info from a JSON file
with open('rois.json', 'r') as json_file:
    region_config = json.load(json_file)

for roi, info in region_config.items():
    x_len = info["xLen"]
    y_len = info["yLen"]
    resolution = info['defaultBin']
    projString = info['projection']

    # Calculate coordinates for the sector polygon
    x_min, x_max = -x_len / 2, x_len / 2
    y_min, y_max = -y_len / 2, y_len / 2

    area_extent = (x_min, y_min, x_max, y_max)
    area_id = regionName
    description = regionName
    proj_id = regionName
    projection = prj4_string_to_dict(projString)
    width = x_len
    height = y_len
    area_def = AreaDefinition(area_id, description, proj_id, projection, width, height, area_extent)
    print(area_def)
    result = kd_tree.resample_nearest(grid_def, da.values, area_def, radius_of_influence=50000, epsilon=0.5)

Error:

Exception has occurred: _ArrayMemoryError       (note: full exception trace is shown but execution is paused at: <module>)
Unable to allocate 33.8 TiB for an array with shape (37158912000000,) and data type bool
  File "/home/synrad/sw/mlde-rts/shapefile_extract_slice.py", line 74, in <module> (Current frame)
    result = kd_tree.resample_nearest(grid_def,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 33.8 TiB for an array with shape (37158912000000,) and data type bool

Why is it pulling so much memory? Can I define the source grid by just using 1D coords? Is there a way to do that? The shape of the source data is: 741, 1400 and the target region is: extent: (-4032000.0, -2304000.0, 4032000.0, 2304000.0) in meters.

[djhoese]

Side suggestion: You can get rid of the projection string -> dict conversion. AreaDefinitions understand anything pyproj CRS objects understand so WKT, EPSG, PROJ.4 strings, PROJ.4 dictionaries, or CRS objects themselves should all work.

For your GridDefinition, you don't need to define it as a GridDefinition. If these are equally spaced pixels (which they are) then creating an AreaDefinition with the equivalent CRS (projection) and calculating extents in degrees should be fine. Note the extents are the outer edges of the pixels, not the centers of the pixels. This should at least save you the memory for the creation of the GridDefinition, but...what is the width and height of your target Area? Are you sure on the size of your input data? Can you print out all the shapes involved in this code?

[Marston]

@djhoese Thanks for the quick response and tips. I removed the function and applied your suggestions but got the same error.

updated code:

import os
os.environ["PROJ_LIB"] = "/root/anaconda3/envs/exp/share/proj"
import matplotlib.pyplot as plt
import json
import xarray as xr
import numpy as np
from pyresample.geometry import AreaDefinition, GridDefinition
from pyresample import kd_tree, save_quicklook

# Load the data array from the netCDF file
da = xr.open_dataarray('/root/gfs_orography.nc', engine='netcdf4')
shape = da.shape
area_id = 'globe'
proj_string = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
area_extent = (-180.0, -90.0, 180.0, 90.0)
shape = shape
src_def = AreaDefinition.from_extent(area_id, proj_string, shape, area_extent)

# Load the region projection info from a JSON file
with open('gsw_rts_map_coords/domains.json', 'r') as json_file:
    region_config = json.load(json_file)

for regionName, info in region_config.items():
    # Extract parameters from the definition
    x_len = info["xLen"]
    y_len = info["yLen"]
    resolution = info['defaultBin']
    projString = info['projection']
    # Calculate coordinates for the sector polygon
    x_min, x_max = -x_len / 2, x_len / 2
    y_min, y_max = -y_len / 2, y_len / 2

    area_extent = (x_min, y_min, x_max, y_max)
    area_id = regionName
    description = f'GSWR {regionName}'
    proj_id = regionName
    projection = projString
    width = x_len
    height = y_len
    targ_def = AreaDefinition(area_id, description, proj_id, projection, width, height, area_extent)
    
    result = kd_tree.resample_nearest(src_def, da.values, targ_def, radius_of_influence=50000, epsilon=0.5)

But I think the problem is the source data size:

Area ID: my_area
Description: my_area
Projection ID: my_area
Projection: {'a': '6370997.0', 'ellps': 'sphere', 'lat_0': '50', 'lon_0': '-135', 'over': 'True', 'proj': 'laea', 'type': 'crs', 'units': 'm'}
Number of columns: 8064000
Number of rows: 4608000
Area extent: (-4032000.0, -2304000.0, 4032000.0, 2304000.0)

I see that PyResample is trying to make a grid of the target grid. Can this be avoided?

[djhoese]

What is the my_area you've printed out in the above comment? You said your source data is 721 rows by 1440 columns, but the area you've printed out with the comment "I think the problem is the source data size" has 4608000 rows by 8064000 columns. Your source data and source area should be the same size.

Additionally, that number of rows and columns is HUGE:

In [1]: 8064000 * 4608000
Out[1]: 37158912000000

That's 37 trillion elements. For 32-bit floats (x and y coordinates usually need 64-bit) that's something like 135TB's per coordinate array...right? That target area definition has a resolution of 1 meter per pixel:

In [5]: (4032000.0 - -4032000.0) / 8064000
Out[5]: 1.0

Could you explain again what you're trying to do?

[djhoese]

And no, pyresample needs to compute all the coordinates for these areas so it can:

1. Build a KDTree of all input pixel coordinates
2. Query the KDTree to say "what pixel in the input data is nearest to this output pixel"

[Marston]

Oh man! I forgot to divide the weight and width by the resolution. Now it works like a charm 👍 I was so focused on understanding PyResample that I forgot that part.
Thanks for the extra eyes and tips.

[Marston]

@djhoese

One last question: can save_quicklook() add the labels to meridians and parallels in the plot?

[djhoese]

Glad it was a simple fix. I was really concerned you actually needed that size/resolution.

I honestly didn't know save_quicklook existed. It looks like it is just using cartopy so it should be possible, but I don't see any option for it. You may be able to take the code here and modify it as needed (or make a pull request):

pyresample/pyresample/plot.py

Lines 241 to 269 in 78b19ca

	def _get_quicklook(area_def, data, vmin=None, vmax=None,
	label='Variable (units)', num_meridians=45,
	num_parallels=10, coast_res='110m', cmap='RdBu_r'):
	"""Get default cartopy matplotlib plot."""
	import matplotlib.pyplot as plt
	coast_res, is_cartopy = _translate_coast_resolution_to_cartopy(coast_res)
	if not is_cartopy:
	return _basemap_get_quicklook(
	area_def, data, vmin, vmax, label, num_meridians,
	num_parallels, coast_res=coast_res, cmap=cmap)
	if area_def.shape != data.shape:
	raise ValueError('area_def shape %s does not match data shape %s' %
	(list(area_def.shape), list(data.shape)))
	crs = area_def.to_cartopy_crs()
	ax = plt.axes(projection=crs)
	ax.coastlines(resolution=coast_res)
	ax.set_global()
	if num_meridians or num_parallels:
	_ = _add_gridlines(ax, num_meridians, num_parallels)
	if not (np.ma.isMaskedArray(data) and data.mask.all()):
	col = ax.imshow(data, transform=crs, extent=crs.bounds,
	origin='upper', vmin=vmin, vmax=vmax, cmap=cmap)
	plt.colorbar(col, shrink=0.5, pad=0.05).set_label(label)
	return plt

[Marston]

@djhoese

I think I have solved the grid lines labels but I have to do the PR next week. I have another issue.
It looks like PyResample has issues when the target grid crosses the dateline. This, if memory serves me, is a know issue. Is there a solution to handle projecting across the dateline?

[Marston]

@djhoese

I've gotten the NN to work very well, but ran into a problem with gradient_search for bilinear interpolation. In my case I'm going from Areadefinition to Areadefinition with an input data that is 2D (same as above and which works for the NN example). Here's my function:

from pyresample.geometry import AreaDefinition
from pyresample import kd_tree, gradient
from pyresample import load_area

def gradient_bilinear(source_def, target_def, data):
    resampler = gradient.create_gradient_search_resampler(source_def, target_def)
    regrid = resampler(data)
    return regrid

print(f'Extract parameters from the definition: {regionName}')
targ_def = load_area('/home/rois.yaml', regionName)

# convert to a double
data = da_src.astype('float64')
result = gradient_bilinear(src_def, targ_def,  data.values)

But got the following error:

Exception has occurred: TypeError
'ResampleBlocksGradientSearchResampler' object is not callable
  File "sector_reprojection.py", line 32, in gradient_bilinear
    regrid = resampler(data)
             ^^^^^^^^^^^^^^^

Am I doing this correct? Not much help in the docs on this method.

[Marston]

If I add over=True, it fails for the KD Tree NN method. But as you have seen, it only solves part of the problem.

[mraspaud]

I think I found the problem for gradient resampling: the resampler uses shapely to find the overlap between the source area and destination area. However, shapely doesn't work well with the data line in particular, so the fact the source area is defined between -180 and 180 degrees is a problem. I guess these functions need reworking: https://github.com/pytroll/pyresample/blob/main/pyresample/slicer.py#L152-L188 . I think @ghiggi is working on a fix for this though.

A workaround for this would be to "re-center" the data somehow around the dateline and use not degrees but meters for the area extend units.

[Marston]

@mraspaud Thanks for looking into this. I understand about the dateline, and I can easily keep the source data in a 0:360 longitude grid. But I'm not sure what you mean by the units. How do this factor into the problem?

[mraspaud]

I'm not sure pyresample will handle 0:360 gracefully, as it wants to have -180:180. Hence using plate-carree with meters as a unit should work around the problem, I thin.

[Marston]

Let me see if I am understanding correctly: it the input data in on a plate-carree grid and expressed in meters not degrees, then this resampling with reprojection (LAEA) should work?