Pyresample 1.25.1 create_area_def return wrong lons with the .get_lonlats()

[Elgyii]

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

# pyresample 1.25.1
In [1]: from pyresample import create_area_def

In [2]: a = create_area_def(area_id="nowpap",
   ...:                     projection={'datum': 'WGS84', 'lat_0': '36',
   ...:                     'lon_0': '137', 'proj': 'latlong'},
   ...:                     shape=(1537, 2068),
   ...:                     area_extent=(120, 23, 154, 49))

In [3]: a.get_lonlats()
Out[3]:
(array([[-102.992, -102.975, -102.959, ...,  -69.041,  -69.025,  -69.008],
        [-102.992, -102.975, -102.959, ...,  -69.041,  -69.025,  -69.008],
        [-102.992, -102.975, -102.959, ...,  -69.041,  -69.025,  -69.008],
        ...,
        [-102.992, -102.975, -102.959, ...,  -69.041,  -69.025,  -69.008],
        [-102.992, -102.975, -102.959, ...,  -69.041,  -69.025,  -69.008],
        [-102.992, -102.975, -102.959, ...,  -69.041,  -69.025,  -69.008]]),
 array([[48.992, 48.992, 48.992, ..., 48.992, 48.992, 48.992],
        [48.975, 48.975, 48.975, ..., 48.975, 48.975, 48.975],
        [48.958, 48.958, 48.958, ..., 48.958, 48.958, 48.958],
        ...,
        [23.042, 23.042, 23.042, ..., 23.042, 23.042, 23.042],
        [23.025, 23.025, 23.025, ..., 23.025, 23.025, 23.025],
        [23.008, 23.008, 23.008, ..., 23.008, 23.008, 23.008]]))

Problem description

The latest version of pyresample returns lons out of the area extent.
So, the resampling does not end successfully.
projection_x_coords, projection_y_coords return the correct information though.

a.projection_x_coords
Out[8]: array([120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992])

In [9]: a.projection_y_coords
Out[9]: array([48.992, 48.975, 48.958, ..., 23.042, 23.025, 23.008])

Expected Output

# pyresample 1.24.0
In [1]: from pyresample import create_area_def

In [2]: a = create_area_def(area_id="nowpap",
   ...:                     projection={'datum': 'WGS84', 'lat_0': '36',
   ...:                     'lon_0': '137', 'proj': 'latlong'},
   ...:                     shape=(1537, 2068),
   ...:                     area_extent=(120, 23, 154, 49))

In [3]: a.get_lonlats()
Out[3]:
(array([[120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992],
        [120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992],
        [120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992],
        ...,
        [120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992],
        [120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992],
        [120.008, 120.025, 120.041, ..., 153.959, 153.975, 153.992]]),
 array([[48.992, 48.992, 48.992, ..., 48.992, 48.992, 48.992],
        [48.975, 48.975, 48.975, ..., 48.975, 48.975, 48.975],
        [48.958, 48.958, 48.958, ..., 48.958, 48.958, 48.958],
        ...,
        [23.042, 23.042, 23.042, ..., 23.042, 23.042, 23.042],
        [23.025, 23.025, 23.025, ..., 23.025, 23.025, 23.025],
        [23.008, 23.008, 23.008, ..., 23.008, 23.008, 23.008]]))

Actual Result, Traceback if applicable

Versions of Python, package at hand and relevant dependencies

pyresample 1.25.1

[djhoese]

It depends what we consider "wrong", but obviously if resampling is no longer working that's a problem. In #447 I actually made an attempt at fixing some assumptions and shortcuts that pyresample has had for a long time. These fixes generally mean more logic is given to the pyproj python library and PROJ C library that pyproj wraps.

I think the main thing you are running into is the removal of a shortcut where pyproj transformations of geographic (lon/lat) datasets were not actually given to pyproj and we assumed that what went in (lon/lat degrees) should come out (lon/lat degrees). The problem with this is it is not accurate when the target projection has a prime meridian or longitude shift. So your case is the exact case I was trying to fix.

Another way to look at that is that pyresample tries to use "standard" lon/lat degrees where 0 longitude is the prime meridian at Greenwich. Everything else is in the project units of the AreaDefinition being used. I can't promise this is consistent with SwathDefinitions since it is something we're trying to be more strict about lately.

Do you have a more complete example of what resampling interfaces you are using? A complete but minimal example that I can use to test and verify things would be great. Especially if it doesn't depend on any data files (random data is fine).

[djhoese]

I should clarify, the get_lonlats is returning the "standard" lon/lats that pyresample expects (think EPSG:4326). The lon/lats you want are actually in the projection space of your AreaDefinition (whose "projection" happens to be a geographic lon/lat projection).

Side note: Be careful defining projections as "latlong" with modern PROJ. PROJ is much more strict about axis ordering and since you put "lat" first it may return latitudes first when pyresample (and you) may expect the longitude value. Specifying longlat should work in more cases.

[Elgyii]

Dear David Hoese
Thanks for the feedback

Do you have a more complete example of what resampling interfaces you are using?

# pyresample==1.25.1
In [1]: from pyresample import create_area_def
   ...: from pyresample import geometry
   ...: from pyresample.kd_tree import resample_nearest
   ...: import numpy as np

In [2]: a = create_area_def(area_id="nowpap",
   ...:                     projection={'datum': 'WGS84', 'lat_0': '36',
   ...:                                 'lon_0': '137', 'proj': 'latlong'},
   ...:                     shape=(1537, 2068),
   ...:                     area_extent=(120, 23, 154, 49))
   ...:
   ...: x = np.linspace(120.0, 154.0, 3000)
   ...: y = np.linspace(23, 49, 2000)
   ...: x, y = np.meshgrid(x, y)
   ...: src_geo = geometry.SwathDefinition(x, y)
   ...:
   ...: s = np.sin(np.linspace(0, 1, x.size)).reshape(x.shape)
   ...: out = resample_nearest(
   ...:     source_geo_def=src_geo,
   ...:     data=s,
   ...:     target_geo_def=a,
   ...:     fill_value=None,
   ...:     radius_of_influence=30000
   ...: )
   ...: out[~out.mask]
Out[2]:
masked_array(data=[],
             mask=[],
       fill_value=1e+20,
            dtype=float64)

# pyresample==1.24.0
In [1]: from pyresample import create_area_def
   ...: from pyresample import geometry
   ...: from pyresample.kd_tree import resample_nearest
   ...: import numpy as np

In [2]: a = create_area_def(area_id="nowpap",
   ...:                     projection={'datum': 'WGS84', 'lat_0': '36',
   ...:                                 'lon_0': '137', 'proj': 'latlong'},
   ...:                     shape=(1537, 2068),
   ...:                     area_extent=(120, 23, 154, 49))
   ...:
   ...: x = np.linspace(120.0, 154.0, 3000)
   ...: y = np.linspace(23, 49, 2000)
   ...: x, y = np.meshgrid(x, y)
   ...: src_geo = geometry.SwathDefinition(x, y)
   ...:
   ...: s = np.sin(np.linspace(0, 1, x.size)).reshape(x.shape)
   ...: out = resample_nearest(
   ...:     source_geo_def=src_geo,
   ...:     data=s,
   ...:     target_geo_def=a,
   ...:     fill_value=None,
   ...:     radius_of_influence=30000
   ...: )

In [3]: out[~out.mask]
Out[3]:
masked_array(
  data=[[[8.4093046e-01, 8.4093052e-01, 8.4093070e-01, ...,
          8.4120035e-01, 8.4120053e-01, 8.4120065e-01],
         [8.4065974e-01, 8.4065986e-01, 8.4066004e-01, ...,
          8.4092993e-01, 8.4093010e-01, 8.4093016e-01],
         [8.4038889e-01, 8.4038895e-01, 8.4038913e-01, ...,
          8.4065920e-01, 8.4065938e-01, 8.4065950e-01],
         ...,
         [1.5001664e-03, 1.5003331e-03, 1.5006664e-03, ...,
          1.9991656e-03, 1.9994990e-03, 1.9996658e-03],
         [1.0001666e-03, 1.0003333e-03, 1.0006666e-03, ...,
          1.4991664e-03, 1.4994997e-03, 1.4996664e-03],
         [5.0016673e-04, 5.0033338e-04, 5.0066673e-04, ...,
          9.9916663e-04, 9.9950004e-04, 9.9966663e-04]]],
  mask=False,
  fill_value=3.4028235e+38,
  dtype=float32)

Side note: Be careful defining projections as "latlong" with modern PROJ. PROJ is much more strict about axis ordering and since you put "lat" first it may return latitudes first when pyresample (and you) may expect the longitude value. Specifying longlat should work in more cases.

Thanks for this side note.

Maybe I should have said "inconsistent behaviour between different different pyresample versions"?

[djhoese]

Thanks for the example code, it really helps. There are some small mistakes here that I think I'm starting to understand.

If I take your example I get the same results as is. If I change the projection to longlat, same results (so that's not the problem). If I subtract 137 from your two X extents so they are:

area_extent=(120 - 137, 23, 154 - 137, 49)

then you get the same results as earlier pyresample versions. Note I changed your final statement to print(out[~out.mask].size) to compare, not the actual values. This makes sense in that by specifying lon_0 you are defining the origin of your projection space (0, 0) to be at EPSG:4326 (137E, 0N). Note lat_0 has no effect on a latlong projection. Apparently this "projection" has to be at the equator.

We can further analyze this by using pyproj in ipython:

In [1]: from pyproj import Transformer, CRS

In [2]: a = CRS.from_user_input(4326)

In [3]: b = CRS.from_user_input({'datum': 'WGS84', 'lat_0': '36', 'lon_0': '137', 'proj': 'latlong'})

In [4]: c = CRS.from_user_input({'datum': 'WGS84', 'lat_0': '52', 'lon_0': '137', 'proj': 'latlong'})

In [5]: Transformer.from_crs(a, b).transform(0, 0)
Out[5]: (-137.0, 0.0)

In [6]: Transformer.from_crs(a, c).transform(0, 0)
Out[6]: (-137.0, 0.0)

In [7]: Transformer.from_crs(a, b).transform(0, 15)
Out[7]: (-122.0, 0.0)

In [8]: Transformer.from_crs(a, c).transform(0, 15)
Out[8]: (-122.0, 0.0)

In [9]: Transformer.from_crs(a, b, always_xy=True).transform(0, 15)
Out[9]: (-137.0, 14.999999999999998)

In [10]: Transformer.from_crs(a, c, always_xy=True).transform(0, 15)
Out[10]: (-137.0, 14.999999999999998)

So the CRSes are lon/lat (EPSG:4326), your desired CRS, and your desired CRS but with a different lat_0. In the first two calls we see that even though b and c are different in the definition the transform of (0, 0) EPSG:4326 to the desired projections yields the same result.

In the next two calls we see that they are producing the same results still...but we changed the latitude input, right? We should have expected the "X" value to stay the same and the latitude/Y go up by 15. Well turns out EPSG:4326 is defined with axes (lat, lon). Apparently I was wrong with proj=latlong:

In [12]: a.axis_info
Out[12]:
[Axis(name=Geodetic latitude, abbrev=Lat, direction=north, unit_auth_code=EPSG, unit_code=9122, unit_name=degree),
 Axis(name=Geodetic longitude, abbrev=Lon, direction=east, unit_auth_code=EPSG, unit_code=9122, unit_name=degree)]

In [13]: b.axis_info
Out[13]:
[Axis(name=Longitude, abbrev=lon, direction=east, unit_auth_code=EPSG, unit_code=9122, unit_name=degree),
 Axis(name=Latitude, abbrev=lat, direction=north, unit_auth_code=EPSG, unit_code=9122, unit_name=degree)]

In [14]: c.axis_info
Out[14]:
[Axis(name=Longitude, abbrev=lon, direction=east, unit_auth_code=EPSG, unit_code=9122, unit_name=degree),
 Axis(name=Latitude, abbrev=lat, direction=north, unit_auth_code=EPSG, unit_code=9122, unit_name=degree)]

By adding always_xy=True we can tell pyproj, ignore CRS axis order, just assume (x, y) meaning (lon, lat). Either way, we see that the lat_0 parameter isn't having any effect on the projection.

Bottom line

Your extents in your original AreaDefinition are incorrect. They should be in the projection space of your area. In your case that means the origin is at EPSG:4326 (137, 0). So if you want EPSG:4326 longitude 130 that would be -7 in your AreaDefinition.

...and remove lat_0 from your projection definition.

What do you think? Did I make any sense here?

[Elgyii]

Definitely it makes a lot of sense.
It helps clarifying a lot of issues. Many thanks

However, by saying “Your extents in your original AreaDefinition are incorrect” we are back to your first saying of WHAT YOU CONSIDER INCORRECT.

If I use projected coordinates, then pyresample 1.24.0 behaves likes 1.25.1, which means print(out[~out.mask].size) is zero.

In the language of projections, definitely the latest version is right.

But the issue I had was this change in pyresample behavior.
Imagine this update to latest pyresample, no more results coming!

When I first encountered it, I went to pyresample documentation page to try to understand it, but did not find any discussion about this behavior change.

From your example, I understand that if I set lon_0=0 in the proj dict then my problems are resolved. Since lat_0 is always 0.

But I think it would be nice to have somewhere in the docs page a discussion about the new behavior which is different from the past versions. Or if I missed it, I would appreciate it if you could direct me to that information.
What do you think?
Thanks for helping me clear this up.

[djhoese]

You make very good points. The pyresample is extremely lacking and in my opinion barely useful for modern pyresample. It is one of my goals once I eventually get some time to rewrite the documentation (see #434).

The best I can give you is the CHANGELOG that is not linked to in the web site (but should be) that mentions "fixing" lon/lat conversion when prime meridian is not 0: https://github.com/pytroll/pyresample/blob/main/CHANGELOG.md. Agreed that this is not as public as it should be. Dealing with lon/lat transformations that don't have a 0 prime meridian is a pretty new concept to me...that's my best excuse for the lack of information on it.

Sorry for all the confusion.

[Elgyii]

Thanks a lot David

I really appretiate the pyresample package.
It makes my work easier.
I give trainings on ocean colour data analysis and I use pyresample a lot.

If there is anything I can contribute back, I would be happy to do it.

[Elgyii]

Dear David,

I think the longlat "projection" is quite problematic?

I run into another issue. It might be covered somewhere.
In that case I apologise for the repetition.

When I call create_area_def

1. without providing shape information, I got a bunch of errors.
   The failure points to the fact that resolution in metres (say 250) result in values out of bounds.

However,
2. when I pass the resolution in degrees (say 0.00225), all goes when but the this resolution is reprorted back by pyresample in metres. This behaviour is never there when using other projections, say, laea projection for instance.

Min example:

extent= (-10.578110399046421, 27.046055, 12.60629960095357, 46.45152)
proj= {'datum': 'WGS84', 'lat_0': 36.96185259015367, 'lon_0': 133.78974039904642, 'proj': 'latlong'}
AreaID= 'nowpap_region'
a = create_area_def(projection=proj,
                    area_extent=extent,
                    resolution=0.00225,
                    units='degrees',
                    area_id=AreaID)

WARNING:root:shape found from radius and resolution does not contain only integers: (8624.651111124404, 10304.182222205558) Rounding shape to (8625, 10305) and resolution from (0.002250000000003638, 0.0022499999999965326) meters to (0.0022498214459000477, 0.0022499089855072467) meters

The reason I am writing this is because I work with swath images and most of the time the resolution is reported in metric system. But if I want to create my area in degrees, then I have to estimate the resolution in degrees so that pyresample does not fail. Alternatively, I should suply the shape and all goes well. But the shape also depends on the resolution.

I thought there should a better way. Appreciate your advice on this.

[djhoese]

The failure points to the fact that resolution in metres (say 250) result in values out of bounds.

What errors/failures?

This behaviour is never there when using other projections, say, laea projection for instance.

Most other projections are in metered space/units. The create_area_def utility function defaults to using whatever units your target projection is in. So if you have a latlong projection then it assumes your extents and resolution parameters are in degrees. If you want them to be in meters then you can pass units="meters" but both extents and resolution would have to be in meters in that case.

The only way currently to express parameters in different units is to wrap them in an xarray DataArray object (I regret asking the contributor who added this function to do it this way, but here we are). So you could leave units="degrees" but then do resolution=xr.DataArray(250, attrs={"units": "meters"}).

[Elgyii]

What errors/failures?

I think is no longer relevant to post it since you explain that if you have a latlong projection then it assumes your extents and resolution parameters are in degrees

So you could leave units="degrees" but then do resolution=xr.DataArray(250, attrs={"units": "meters"})

extent= (-10.578110399046421, 27.046055, 12.60629960095357, 46.45152)
proj= {'datum': 'WGS84', 'lat_0': 36.96185259015367, 'lon_0': 133.78974039904642, 'proj': 'latlong'}
AreaID= 'nowpap_region'
a = create_area_def(projection=proj,
                    area_extent=extent,
                    resolution=xr.DataArray(250, attrs={"units": "meters"}),
                    units='degrees',
                    area_id=AreaID)

ValueError                                Traceback (most recent call last)
Cell In [17], line 4
      2 proj= {'datum': 'WGS84', 'lat_0': 36.96185259015367, 'lon_0': 133.78974039904642, 'proj': 'latlong'}
      3 AreaID= 'nowpap_region'
----> 4 a = create_area_def(projection=proj,
      5                     area_extent=extent,
      6                     resolution=xr.DataArray(250, attrs={"units": "meters"}),
      7                     units='degrees',
      8                     area_id=AreaID)
File lib\site-packages\pyresample\area_config.py:516, in create_area_def(area_id, projection, width, height, area_extent, shape, upper_left_extent, center, resolution, radius, units, **kwargs)

513 # Fills in missing information to attempt to create an area definition.

514 if area_extent is None or shape is None:

515     area_extent, shape, resolution = 

--> 516         _extrapolate_information(area_extent, shape, center, radius,

517                                  resolution, upper_left_extent, units,

518                                  p, crs)

519 return _make_area(area_id, description, proj_id, projection, shape,

520                   area_extent, resolution=resolution, **kwargs)
File lib\site-packages\pyresample\area_config.py:777, in _extrapolate_information(area_extent, shape, center, radius, resolution, upper_left_extent, units, p, crs)

775     radius = _convert_units(radius, 'radius', units, p, crs, center=center)

776 # Convert resolution to meters if given as an angle. If center is not found, an exception is raised.

--> 777 resolution = _convert_units(resolution, 'resolution', units, p, crs, center=center)

778 # Inputs unaffected by data below: area_extent is not an input. However, output is used below.

779 if radius is not None and resolution is not None:

780     # Function 2-A
File lib\site-packages\pyresample\area_config.py:656, in _convert_units(var, name, units, p, crs, inverse, center)

654     var = tuple(var.data.tolist())

655 if crs.is_geographic and not ('deg' == units or 'degrees' == units):

--> 656     raise ValueError('latlon/latlong projection cannot take {0} as units: {1}'.format(units, name))

657 # Check if units are an angle.

658 is_angle = ('deg' == units or 'degrees' == units)
ValueError: latlon/latlong projection cannot take meters as units: resolution

[djhoese]

Ugh sorry, it is possible to do the opposite (a metered projection using degree extents). The only way to map some meters to a value in degrees semi-accurately is to have some projection associated with it. We had decided not to "just pick one" for this operation in create_area_def due to the confusion and possible inaccuracy it could produce. I'm sorry I didn't realize this was your use case earlier.

Basically if you did "250m" at the equator in a mercator projection you'd get one resolution, if you did it at latitude 55 you'd get another, if you did it at the equator of a polar-stereographic projection you'd get yet another result. Having pyresample be smart about this would be only an estimate and ugly/difficult to implement.

Another thought, do you need a lon/lat projection or would an equirectangular projection (+proj=eqc) work? It is the same distortion, just in meters.

[Elgyii]

Thank you!
I totally understand the limitations.

I think I will stick to using the shape of the target area.
Appreciated!