Add pixel interpretation and casting logic (AREA_OR_POINT raster metadata) and add new package-level config defaults

[rhugonnet]

Summary

This PR adds support for pixel interpretation (= raster metadata "AREA_OR_POINT") across all functionalities of GeoUtils.
All functions dealing with Raster-Raster, Raster-Point, or coordinate manipulation now have a logic linked to pixel interpretation:

• Raster-Raster functions (arithmetic, use as match for reproject) raise a warning to the user if the two interpretations differs and proposes the user to either redefine one interpretation with set_area_or_point or ignore the warning by changing a global warning configuration (details below),
• Coordinate manipulation functions (coords, xy2ij, ij2xy) automatically shift by half a pixel according to the interpretation, can also be turned off by a global configuration,
• Raster-Point functions utilize coordinate functions and thus also shift according to the interpretation.

For functions that could perform a shift (e.g., interp_points, to_pointcloud, coords), shifting with pixel interpretation can be specified with a boolean argument shift_area_or_pixel.
Additionally, there is a new global configuration dictionary geoutils.config which can be used to change the default behaviour at the package scale. For instance, geoutils.config["shift_area_or_point"]=False removes the shifting according to pixel interpretation everywhere, and geoutils.config["warn_area_or_point"]=False removes the warnings associated with pixel interpretation everywhere also.

What is pixel interpretation ?

This is a confusing topic, see references linked in #503 for more info.
In short, the metadata "AREA_OR_POINT" determines the pixel interpretation:

• "Area" (the most common): the pixel value represents something derived across the entire cell (e.g. an average), in which case its coordinate is tied to the upper left corner which is the raster default,
• "Point" (less common, often used in DEMs): the pixel value represents the pixel center, where the coordinate is tied.

GDAL has some internal logic introduced in 2010 to deal with this attribute but it only applies to tie points (GCP): it shifts them by half a pixel when the interpretation is "Point".

Rasterio derives coordinates independently of this metadata, and thus we have to account for it.

Detailed changes

This PR:

• Adds an area_or_point property attribute to Raster, that can be set using set_area_or_point() (optionally shifts the raster) or using area_or_point=... with default parameters,
• Adds specific test for area_or_point,
• Adds a _config.py module with a config.ini file to set global package parameters, relying the configparser package (included as part of core Python), and a related test_config.py test module,
• Adds two package-wide parameters "shift_area_or_point" and "warn_area_or_point", both defaulting to True, which can be changed by the user as simply as geoutils.config[param] = new_val,
• Makes ij2xy() and coords() understand pixel interpretation (only xy2ij() had an implementation until now), now the ij2xy() and xy2ij() are fully reversible for any pixel interpretation and coords() are derived from ij2xy() for full consistency. The offset argument has become a force_offset argument (because coordinates are always defined as upper-left or center depending on interpretation for raster, we should not leave the choice to the user for an offset value except for a "forced" call),
• Updates interp_points() to maintain behaviour consistent with coordinates functions,
• Adds many tests and adapts existing tests for all coordinates functions,
• Adds a warning in reproject() and crop() if a match-reference raster is used that has a different pixel interpretation than the one being reprojected/cropped,
• Adds area_or_point and tags arguments to from_array() for passing any pixel interpretation or tag metadata to a raster defined from array, and modifies copy() accordingly,
• Changes all calls of from_array() that only changed data= into calls to copy(new_array=) (completely equivalent and no need to pass on the area_or_point),
• Adds specific tests to check that self.ds.tag()['AREA_OR_POINT'] gets updated into 'Area' when cropping (updating) #127 is resolved,
• Merges Raster._overloading_check and _check_cast_array_raster into a single _cast_numeric_array_raster function,
• Opens issues Nodata should be cast depending on operation type and not just input #517 and User-input type checks of nodata bypassed in from_array() and other dependent functions #516 in relation to nodata casting.

Resolves #508
Resolves #500
Resolves #503
Resolves #127

TODO

• Merge casting functions following Merge _check_cast_raster_array() and _overloading_check() into a consistent implementation #508,
• Add shift_area_or_point transformation to all coordinates function still missing it (coords and ij2xy),
• Add warning on area_or_point when using a match-reference raster to crop/reproject/proximity another raster?
• Modify shift_area_or_point to be None by default and call config["shift_area_or_point"],
• Adjust interp_points functionality,
• Add option to convert interpretation with shift in set_area_or_point(),
• Remove forced "Area" interpretation for outputs of crop(),
• Add tests on direction of shift for both coordinates and transform.

Issues for another PR

Making a documentation page on pixel interpretation, opened #523

The issues #516 and #517 need to be addressed for consistent output, and will require additional modifications of the merged casting function _cast_numeric_array_raster.

[rhugonnet]

@adehecq @erikmannerfelt The new config could also be used to set default matplotlib interpolation, or number of threads used during reprojecting 😜
Which we discussed in GlacioHack/xdem#383 and #438

[adehecq]

A welcome clarification of the pixel interpretation ! 👏 Thanks !
It's still not fully clear in my head why the signs of the half pixels shifts are sometimes different, but maybe you'll be able to help answer.
I'm also curious how often this pixel interpretation will matter. I rarely encountered a case when I needed to care about it.
Oh and maybe this would require adding a small section about this in the documentation?

[rhugonnet]

Thanks for the review!
I clarified why the half-pixel shift sign differs between the two functions and added a lot of tests in 4cf8899. All working as intended! 🥳 (based on the references in #503)
And opened issue #523 to add a page in the documentation soon.

And for the use: Basically, you should never have to touch the default parameters, except if you get a warning 😄