Implement distance_to_anomaly (Closes #209) #324

nialov · 2024-02-13T15:58:04Z

Performance issues are a problem which can be solved by using the alternative functionality that uses gdal_proximity function instead of existing distance_computation implementation. The GDAL implementation is orders of magnitude faster.

Draft for now until I clean up the implementation.

@tomironkko made an initial draft of this function which might have been better than the current shapely distance calculation that uses the existing distance_computation one but more than likely much worse than the gdal_proximity method.

nialov · 2024-02-14T10:04:07Z

Ready for review. The main function, distance_to_anomaly, works on all platforms while distance_to_anomaly_gdal only works on Linux and MacOS. Based on initial tests I thought the GDAL utility functions would work on all platforms but found out later they are inconsistent on Windows. The main problem on Windows is the lack of osgeo_utils module in the conda environment.

Anyway, thought it best to keep the GDAL method for later reference. I just marked the corresponding test as an xfail when run on Windows. It is orders of magnitude faster and at least somewhat more memory-efficient (filtering to anomalous values is still done in-memory).

If distance_computation, which the main function relies on, requires optimization later, you should check out the original branch by @tomironkko here: master...209-add-distance-to-anomaly. The "raster" distance math approach is more efficient but it can obviously only handle raster data while the current distance_computation handles vector data as well. You could rasterize first and then use raster math. Alternatively the missing osgeo_utils could be investigated.

nmaarnio

Hey, looks really good overall! I wrote a few comments that mainly regard docs and checks. I think we can keep the GDAL version.

nmaarnio · 2024-02-15T10:07:55Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+        anomaly_raster: The raster in which the distances
+            to the nearest anomalous value are determined.


In the parameters, anomaly_raster_profile and anomaly_raster_data are separate. I suppose you just forgot to update the docs here?

Yep, updating!

nmaarnio · 2024-02-15T10:09:42Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+    raster_profile: Union[profiles.Profile, dict],
+    anomaly_raster_profile: Union[profiles.Profile, dict],
+    anomaly_raster_data: np.ndarray,


A question about these: Is raster_profile the target/output profile? For me, it was not (and maybe isn't still) completely clear why we need to specify two raster profiles. In what cases raster_profile != anomaly_raster_profile? If we keep both, the naming/docs could be clarified.

Also, our raster processing functions have taken in Rasterio's DatasetReaders so far and this implementation conflicts that convention. However, I have to say I prefer this way in general, so maybe this does not need to be changed and we can considering refactoring the other tools.

With the initial (non-gdal) approach it was possible to calculate the distances for a raster that is not equivalent in profile to the raster in which the anomalies are determined from. So anomaly raster could have e.g. different resolution to the raster in which the distances were calculated to. For simplification sense I have now removed this option, I suppose it is probably not a common use case. Easier to maintain without it 👍

nmaarnio · 2024-02-15T10:11:30Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+    threshold_criteria_value: THRESHOLD_CRITERIA_VALUE_TYPE,
+    threshold_criteria: THRESHOLD_CRITERIA_TYPE,


These input types do not transfer to the docs well. I would rather just put the Union/Literal types directly to the function signature for better user experience and documentation even if I understand they are used repeatedly in this module.

Yeah I keep trying to minimize code related to type annotations but Python does not make it easy 😢 . I will update to the explicit versions.

nmaarnio · 2024-02-15T10:12:12Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+    anomaly_raster_data: np.ndarray,
+    threshold_criteria_value: THRESHOLD_CRITERIA_VALUE_TYPE,
+    threshold_criteria: THRESHOLD_CRITERIA_TYPE,
+) -> np.ndarray:


We have had the convention of always returning raster profile/meta in raster processing tools. I would follow this schema in this tool too.

nmaarnio · 2024-02-15T10:17:24Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+    raster_width = raster_profile.get("width")
+    raster_height = raster_profile.get("height")
+
+    if not isinstance(raster_width, int) or not isinstance(raster_height, int):
+        raise InvalidParameterValueException(
+            f"Expected raster_profile to contain integer width and height. {raster_profile}"
+        )
+
+    raster_transform = raster_profile.get("transform")
+
+    if not isinstance(raster_transform, transform.Affine):
+        raise InvalidParameterValueException(
+            f"Expected raster_profile to contain an affine transformation. {raster_profile}"
+        )


Is it on purpose that we only check raster_profile but not anomaly_raster_profile?

While reviewing this I got an idea: should we add a check called check_raster_meta or check_raster_profile that allows selecting (or deselecting) profile/meta elements to check? This check could essentially do what you are doing here for the specified elements and raise the exceptions too. Currently many of our raster processing functions do not check if metadata contain the required elements.

If we add such a check (and maybe even if we didn't), we could also consider a new dedicated exception type such as InvalidRasterMetaException.

em-t · 2024-02-15T16:31:13Z

I did some testing with the non-gdal version using various raster files and using a randomly generated anomaly profile. Everything has looked good so far 👍 I'm wondering, though, if the function should handle nodata values in some more sophisticated way, or is it ok to assume that the user has replaced nodata values with nan for the input raster?

nialov · 2024-02-23T13:45:15Z

I did some testing with the non-gdal version using various raster files and using a randomly generated anomaly profile. Everything has looked good so far 👍 I'm wondering, though, if the function should handle nodata values in some more sophisticated way, or is it ok to assume that the user has replaced nodata values with nan for the input raster?

Yes, you are correct that if np.nan is not the true nodata value, then the nodata values will be left among the possible anomaly values. I have personally tried to advocate for a centralized approach to handling of nodata but if there are no current plans for such approach, I can implement nodata handling in this function as well (@nmaarnio)?

nialov · 2024-02-23T13:47:21Z

Ready for review, let me know about the nodata handling!

nmaarnio · 2024-02-29T10:42:38Z

Let's implement nodata handling in this tool now. However, if you have time and interest to suggest/implement a centralized approach for nodata handling across the toolkit in the upcoming weeks/months, please go ahead.

Performance issues are a problem which can be solved by using the alternative functionality that uses gdal_proximity function instead of existing distance_computation implementation. The GDAL implementation is orders of magnitude faster.

@nmaarnio

Based on review by @nmaarnio and @em-t, updated faulty documentation, simplified inputs and outputs and made a separate check for raster profile values.

Updated tests according to previous changes and added some new ones to test raster profile check.

Uses nodata value from raster profile and falls back to using np.nan if profile does not contain it.

nialov · 2024-03-07T12:40:14Z

Nodata handling implemented and ready for review again. Nodata only affects the definition of anomalous values in _fits_criteria. The output will never have nodata values in it.

nmaarnio

I noticed one potential error in the code, see my comment below. Also, I think it could be useful to add a check that threshold_criteria_value and threshold_criteria match (to check if threshold_criteria_value is a tuple when "in_between" or "outside" is used – not sure if numeric checks for the threshold_criteria_value should be included). I couldn't test this locally and check the notebook, but reading the code looks otherwise good to me!

nmaarnio · 2024-03-08T09:33:37Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+        "in_between": lambda anomaly_raster_data: np.where(
+            np.logical_and(anomaly_raster_data > threshold_criteria[0], anomaly_raster_data < threshold_criteria[1])
+        ),
+        "outside": lambda anomaly_raster_data: np.where(
+            np.logical_or(anomaly_raster_data < threshold_criteria[0], anomaly_raster_data > threshold_criteria[1])


It looks like threshold_criteria here should bethreshold_criteria_value.

Thanks, fixed and added tests for each criteria!

nialov · 2024-03-11T06:58:12Z

Thanks for the review, @nmaarnio! Ready for review again 👍

nmaarnio

LGTM now!

em-t

Sorry for leaving my review hanging for so many days. I did some more testing and it looks good to me! I left two minor suggestions, but I'm also ok with merging now. Nice job 👍

em-t · 2024-03-14T11:22:08Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+    threshold_criteria_value: Union[Tuple[Number, Number], Number],
+    threshold_criteria: Literal["lower", "higher", "in_between", "outside"],
+    anomaly_raster_data: np.ndarray,
+    nodata_value: Optional[Number],


Should this default to Optional[Number] = None?

(The value of the Optional type-hinted parameter doesn't seem to get set to None by default if omitted. The type hint just allows for the type to be either a number or NoneType.)

Because the function is a private one I think it is best to make sure all inputs are explicitly passed when using the function so no functionality is missed. Using Optional does not mean you need to set the default value, it just means that input is allowed to be passed as None.

em-t · 2024-03-14T11:58:00Z

eis_toolkit/raster_processing/distance_to_anomaly.py

+        for row in rows
+    ]
+    all_points = list(chain(*all_points_by_rows))
+    all_points_gdf = gpd.GeoDataFrame(geometry=all_points, crs=anomaly_raster_profile["crs"])


(Just a suggestion, feel free to ignore this if you think best.)

If all_points ends up empty due to nothing matching the criteria, this function could raise an exception with an informative message. The distance_computation does raise an exception when checking the shape, but the message is sort of generic, and it might be good to let the user know why it's happened in this case.

Thanks, added a check that the criteria matches at least some of the data!

If the criteria does not match any of the data, an exception is raised to notify the user.

nialov · 2024-03-15T07:04:11Z

Ready for review/merging 👍 @nmaarnio, @em-t

em-t · 2024-03-15T08:24:43Z

Looks good!

nmaarnio · 2024-03-15T09:18:23Z

Merging!

nialov force-pushed the 209-add-distance-to-anomaly-v2 branch from 847eb53 to c789da9 Compare February 14, 2024 09:47

nialov marked this pull request as ready for review February 14, 2024 10:04

nialov requested review from nmaarnio, msmiyels and em-t February 14, 2024 10:04

nmaarnio requested changes Feb 15, 2024

View reviewed changes

nialov added 6 commits March 7, 2024 13:35

feat: implement distance_to_anomaly

bbf437d

Performance issues are a problem which can be solved by using the alternative functionality that uses gdal_proximity function instead of existing distance_computation implementation. The GDAL implementation is orders of magnitude faster.

docs: fix names

5b0057c

refactor: use Number as type

1dc9c86

fix(distance_to_anomaly): apply review suggestions

7775bfb

Based on review by @nmaarnio and @em-t, updated faulty documentation, simplified inputs and outputs and made a separate check for raster profile values.

test(distance_to_anomaly): update tests

abd858e

Updated tests according to previous changes and added some new ones to test raster profile check.

docs(distance_to_anomaly): update notebook

5d25fdb

nialov force-pushed the 209-add-distance-to-anomaly-v2 branch from 334d523 to 5d25fdb Compare March 7, 2024 11:35

nialov added 4 commits March 7, 2024 14:21

feat(distance_to_anomaly): handle nodata

ead23e7

Uses nodata value from raster profile and falls back to using np.nan if profile does not contain it.

test(distance_to_anomaly): test with nodata

b706b25

fix(miscellaneous): fix check logic

cc19baa

docs(distance_to_anomaly): add function input

9b24bad

nmaarnio requested changes Mar 8, 2024

View reviewed changes

nialov added 5 commits March 8, 2024 17:53

feat(distance_to_anomaly): check threshold inputs

ed187cb

fix(distance_to_anomaly): use correct variable

abcc7ba

test(distance_to_anomaly): test checks and all criteria

819d5d5

docs(distance_to_anomaly): add all criteria examples

cbcede4

style(distance_to_anomaly): clean code and docs

53e3edd

nmaarnio approved these changes Mar 11, 2024

View reviewed changes

em-t approved these changes Mar 14, 2024

View reviewed changes

nialov added 2 commits March 15, 2024 08:58

feat(distance_to_anomaly): check threshold data match

74cc1da

If the criteria does not match any of the data, an exception is raised to notify the user.

test(distance_to_anomaly): test criteria data check

fea0258

nmaarnio approved these changes Mar 15, 2024

View reviewed changes

nmaarnio merged commit e01bf58 into GispoCoding:master Mar 15, 2024
6 checks passed

nialov mentioned this pull request Jun 10, 2024

Optimize Distance to anomaly #409

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement distance_to_anomaly (Closes #209) #324

Implement distance_to_anomaly (Closes #209) #324

nialov commented Feb 13, 2024

nialov commented Feb 14, 2024

nmaarnio left a comment

nmaarnio Feb 15, 2024

nialov Feb 23, 2024

nmaarnio Feb 15, 2024

nialov Feb 23, 2024

nmaarnio Feb 15, 2024

nialov Feb 23, 2024

nmaarnio Feb 15, 2024

nialov Feb 23, 2024

nmaarnio Feb 15, 2024

em-t commented Feb 15, 2024

nialov commented Feb 23, 2024

nialov commented Feb 23, 2024

nmaarnio commented Feb 29, 2024

nialov commented Mar 7, 2024

nmaarnio left a comment

nmaarnio Mar 8, 2024

nialov Mar 8, 2024

nialov commented Mar 11, 2024

nmaarnio left a comment

em-t left a comment

em-t Mar 14, 2024

nialov Mar 15, 2024

em-t Mar 14, 2024

nialov Mar 15, 2024

nialov commented Mar 15, 2024

em-t commented Mar 15, 2024

nmaarnio commented Mar 15, 2024

		anomaly_raster: The raster in which the distances
		to the nearest anomalous value are determined.

		threshold_criteria_value: THRESHOLD_CRITERIA_VALUE_TYPE,
		threshold_criteria: THRESHOLD_CRITERIA_TYPE,

Implement distance_to_anomaly (Closes #209) #324

Implement distance_to_anomaly (Closes #209) #324

Conversation

nialov commented Feb 13, 2024

nialov commented Feb 14, 2024

nmaarnio left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

em-t commented Feb 15, 2024

nialov commented Feb 23, 2024

nialov commented Feb 23, 2024

nmaarnio commented Feb 29, 2024

nialov commented Mar 7, 2024

nmaarnio left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nialov commented Mar 11, 2024

nmaarnio left a comment

Choose a reason for hiding this comment

em-t left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nialov commented Mar 15, 2024

em-t commented Mar 15, 2024

nmaarnio commented Mar 15, 2024