Naive "rounding to even" affect correctness of `get_array_indices_from_lonlat` and `get_area_slices` in some occasions

Problem description

This issue addresses the current inaccurate computation of integer indices from (float) array indices.
This can cause errors in many methods related to indices retrieval (see below) and possible inaccurate results of get_area_slices (see https://github.com/pytroll/pyresample/issues/371)

The fact is that the method round used to get the integer indices from float indices "rounds to the even".

Rounds to even example

round(0.5) # --> 0 <-- lower (even)
round(1.5) # --> 2 <-- upper (odd)
round(2.5) # --> 2 <-- lower (even) 
round(3.5) # --> 4 <-- upper (odd)

To visualize the problem, let's create a classical WGS84 grid at 0.5 resolution of shape (y=360,x=720).
The valid number of indices on x goes from 0 to 359.
The corners are [(-180.0, 90.0), (180.0, 90.0), (180.0, -90.0), (-180.0, -90.0)]

from pyresample.geometry import AreaDefinition
area_def = AreaDefinition.from_epsg(code=4326, resolution=0.5) 
area_def.height
area_def.width
area_def.shape
area_def.outer_boundary_corners  

Now let's try to retrieve the integer indices of the top and bottom left corners

# Integer indices 
area_def.get_array_indices_from_lonlat(-180, 90)  # (0, 0) --> OK         
area_def.get_array_indices_from_lonlat(-180, -90) # ValueError: Point outside area:( 0.000000 360.000000) ---> BUG 

# Float indices 
area_def.get_array_coordinates_from_lonlat(-180, 90)   # (-0.5, -0.5)
area_def.get_array_coordinates_from_lonlat(-180, -90)  # (-0.5, 359.5)

Employing np.round, lead to rounding of 359.5 (even) upper to 360 ... outside the valid integer indices of the AreaDef (0-359)

np.round(area_def.get_array_coordinates_from_lonlat(-180, -90)) # array([ -0., 360.])

Solution

Replace all round() occurence in geometry.py by the below round_to_indices function

def round_to_indices(x): 
    # First include the start border 
    x[x == -0.5] = 0
    # Then round to the lower (instead of even) 
    decimal_values = x - np.floor(x)
    idx_midpoint = np.where(decimal_values == 0.5)[0]
    if len(idx_midpoint) > 0:
        x[idx_midpoint] = x[idx_midpoint] - 0.0000001
    return np.round(x)  

# Rounding now is performed "correctly" 
round_to_indices(area_def.get_array_coordinates_from_lonlat(-180, -90)) # array([  0., 359.])

# Other examples
arr = np.array([0.5, 1.5, 2.5, 3.5, 4.5])
np.round(arr)         # array([0., 2., 2., 4., 4.])  <-- rounding to even 
round_to_indices(arr) # array([0., 1., 2., 3., 4.])

arr = np.array([-0.51, -0.5, 0, 0.5, 1, 358.5, 359, 359.5, 359.1])
np.round(arr) 
round_to_indices(arr) # change behaviour only for X.5 numbers 

Further notes

Addressing the rounding issue, will automatically solve the same problem arising currently in the deprecated methods lonlat2colrow and get_xy_from_lonlat.

area_def.lonlat2colrow(lons=-180, lats=-90) # ValueError: Point outside area:( 0.000000 360.000000) ---> BUG 
area_def.get_xy_from_lonlat(-180, -90)      # ValueError: Point outside area:( 0.000000 360.000000) ---> BUG 

Start in geometry.py by tracing get_array_indices_from_lonlat and get_area_slices, then check the deprecated lonlat2colrow and get_xy_from_lonlat paths mentioned in the issue. Reproduce the WGS84 0.5-resolution corner cases and verify that midpoint rounding stays within valid indices, including the lower-left and lower-right boundaries.