pytroll / pyresample

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

Fix DynamicAreaDefinition not preserving user's requested resolution #523

Closed djhoese closed 1 year ago

djhoese commented 1 year ago

Closes the new issue described in #517 where the user points out that they requested a resolution of 250.0 but the resulting AreaDefinition does not match that and actually changes between cases. This PR fixes this so in the cases where resolution is provided (and not shape) the resulting area should have pixel_size_x/pixel_size_y very very close to this value except for some floating point precision.

codecov[bot] commented 1 year ago

Codecov Report

Merging #523 (9784190) into main (2b94b70) will increase coverage by 0.01%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main     #523      +/-   ##
==========================================
+ Coverage   94.32%   94.33%   +0.01%     
==========================================
  Files          79       79              
  Lines       12944    12982      +38     
==========================================
+ Hits        12209    12247      +38     
  Misses        735      735              
Flag Coverage Δ
unittests 94.33% <100.00%> (+0.01%) :arrow_up:

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files Coverage Δ
pyresample/test/test_gradient.py 99.00% <ø> (-0.01%) :arrow_down:
pyresample/geometry.py 87.53% <100.00%> (+0.02%) :arrow_up:
pyresample/test/test_geometry_legacy.py 97.68% <100.00%> (+0.03%) :arrow_up:

... and 1 file with indirect coverage changes

coveralls commented 1 year ago

Coverage Status

coverage: 93.939% (+0.09%) from 93.846% when pulling 9784190ff115955bb878b7b9f395df9136f71f18 on djhoese:bugfix-dynamic-res-preserve into 2b94b704365507a3ac9b2de9e9a25e5414792b73 on pytroll:main.

djhoese commented 1 year ago

@mraspaud I have some local changes that make the extents nice round numbers (usually integers, no fractional meters), but it results in a lot of the tests failing because shapes end up being off by 1 with the new logic or extents end up being slightly different. I want to run the code by you before I fix all the failing tests. What do you think:

        if shape:
            height, width = shape
            x_resolution = (corners[2] - corners[0]) * 1.0 / (width - 1)
            y_resolution = (corners[3] - corners[1]) * 1.0 / (height - 1)
            area_extent = (corners[0] - x_resolution / 2,
                           corners[1] - y_resolution / 2,
                           corners[2] + x_resolution / 2,
                           corners[3] + y_resolution / 2)
        else:
            x_resolution, y_resolution = resolution
            half_x = x_resolution / 2
            half_y = y_resolution / 2
            # align extents with pixel resolution
            area_extent = (
                math.floor((corners[0] - half_x) / x_resolution) * x_resolution,
                math.floor((corners[1] - half_y) / y_resolution) * y_resolution,
                math.ceil((corners[2] + half_x) / x_resolution) * x_resolution,
                math.ceil((corners[3] + half_y) / y_resolution) * y_resolution,
            )
            width = int(round((area_extent[2] - area_extent[0]) / x_resolution))
            height = int(round((area_extent[3] - area_extent[1]) / y_resolution))

        return area_extent, width, height

The shape block is basically unchanged. The else is the important part. @yukaribbba your feedback is welcome too.

Basically I kept the width/height calculations essentially the same, but instead of adjusting the corners and then doing half-pixel math, I calculation the extents immediately and round them to the nearest pixel position. This has two benefits to the existing code:

  1. I'm not adding pixels only to the right side of the area, but adding space on both sides.
  2. The "round" extents mean that there is very little chance of floating point precision issues adding 0.000000005 to the resulting resolution. Especially in metered projections where we likely have no fractional parts to the extents.
yukaribbba commented 1 year ago

I have no problem with this. Don't forget to exclude decimals.

djhoese commented 1 year ago

Exclude decimals?

yukaribbba commented 1 year ago

Yeah I mean decimal degrees. Wait... no need for that.

djhoese commented 1 year ago

Yeah we should be fine since we're doing everything based on resolution.

mraspaud commented 1 year ago

@djhoese I'm good with that!