dakk / libweatherrouting

A sailing weather routing library
https://dakk.github.io/libweatherrouting
GNU General Public License v3.0
24 stars 11 forks source link

Upwind navigation returns wrong results #7

Closed dakk closed 2 years ago

dakk commented 3 years ago

As reported by @enricofer :

Start at 09:48 del 22/10/2021 from E25.27515,N35.80502 to E24.73033,N35.81287 (E to W)

test02.grb.tar.gz

It returns a wrong result:

image

We need to create a test case to check if libwr behavior is correct

dakk commented 3 years ago

@enricofer I created a test case for an upwind navigation, and the result seems correct; this is a routing from E to W, with a 20knot wind from W:

image

I've found an error in the geojson transform function (lat and lon are inverted in geojson). I add two new tests which print the geojson:

You can start from here to check algorithm behavior .

Reading your code, I can't figure out what's wrong; btw I think this line is wrong, since you are appending destination as last point (which could be against the wind):

https://github.com/enricofer/wind_forecast_routing/blob/master/wind_forecast_routing_algorithm.py#L306

enricofer commented 3 years ago

Thanks for the verify and the code review @dakk. It clearly appears that the library is running as expected. I fear that the issue is in how the plugin read and pass grib samples.

dakk commented 3 years ago

@enricofer put an assert to debug results from getWindAt; using the grib you passed me it's easy to check if values are retrieved correctly since should be a strong W wind

enricofer commented 3 years ago

I'm having trouble checking results for upwind conditions. I'm modified the test suite to check upwind navigation for all directions: https://github.com/enricofer/libweatherrouting/blob/9a876316a0030e2b631f6458a2f2924aab919d44/tests/linearbestisorouter_test.py#L69 Basically, The test case get a track ciclying from a base point to all quadrants, mocking grib with a constant upwind in reverse track direction. For all directions I'm getting a reasonable result but for upwind navigation in north direction I'm getting this not expected result: https://gist.github.com/enricofer/f064c587e13588ef660a4ecf89577d8b These are all test output together in one geojson: https://gist.github.com/enricofer/411843c5d1cd9e73292d84ebf9880358 I don't understand why.

enricofer commented 3 years ago

Furthermore, is it right this result? I included a test on routage speed, I would expect zero for twa=0, as for getSpeed method, but I get different result....

Traceback (most recent call last):
File "/home/travis/build/enricofer/libweatherrouting/tests/polar_test.py", line 36, in test_routage
self.assertEqual(self.polar_obj.getRoutageSpeed(10,math.radians(0)), 0, 0.001)
AssertionError: 4.443057770090072 != 0 : 0.001
dakk commented 2 years ago

Sorry for the delay; I tried your tests and results are wrong as you pointed out. Atm I've no time for a deep debug, but try this: here https://github.com/enricofer/libweatherrouting/blob/9a876316a0030e2b631f6458a2f2924aab919d44/weatherrouting/routers/router.py#L88 instead of using getRoutageSpeed, use getSpeed; your test seems to give a correct result:

{"type": "FeatureCollection", "features": [{"type": "Feature", "id": 1, "geometry": {"type": "Point", "coordinates": [17.0, 34.5]}, "properties": {"strait": "mean"}}, {"type": "Feature", "id": 1, "geometry": {"type": "LineString", "coordinates": [[17.05, 34.5], [18.0, 35]]}, "properties": {"strait": "north"}}, {"type": "Feature", "id": 1, "geometry": {"type": "LineString", "coordinates": [[16.95, 34.5], [16.0, 34]]}, "properties": {"strait": "south"}}, {"type": "Feature", "id": 0, "geometry": {"type": "Point", "coordinates": [17, 34]}, "properties": {"timestamp": "2021-04-02 12:00:00", "attr1": 0, "attr2": 0, "attr3": 0, "attr4": 0}}, {"type": "Feature", "id": 1, "geometry": {"type": "Point", "coordinates": [16.887493858644547, 34.13356039140768]}, "properties": {"timestamp": "2021-04-02 13:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 5.275, "attr4": 325.0}}, {"type": "Feature", "id": 2, "geometry": {"type": "Point", "coordinates": [17.050774167231765, 34.269009846086995]}, "properties": {"timestamp": "2021-04-02 14:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 6.2, "attr4": 45.0}}, {"type": "Feature", "id": 3, "geometry": {"type": "Point", "coordinates": [16.93790920922015, 34.40256384863036]}, "properties": {"timestamp": "2021-04-02 15:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 5.275, "attr4": 325.0}}, {"type": "Feature", "id": 4, "geometry": {"type": "Point", "coordinates": [17.101712917701974, 34.53800624320647]}, "properties": {"timestamp": "2021-04-02 16:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 6.2, "attr4": 45.0}}, {"type": "Feature", "id": 5, "geometry": {"type": "Point", "coordinates": [16.962186228592657, 34.67529960820302]}, "properties": {"timestamp": "2021-04-02 17:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 5.8, "attr4": 320.0}}, {"type": "Feature", "id": 6, "geometry": {"type": "Point", "coordinates": [17.126527752866963, 34.81073481562158]}, "properties": {"timestamp": "2021-04-02 18:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 6.2, "attr4": 45.0}}, {"type": "Feature", "id": 7, "geometry": {"type": "Point", "coordinates": [16.98654060783366, 34.948021207463654]}, "properties": {"timestamp": "2021-04-02 19:00:00", "attr1": 0.0, "attr2": 10.0, "attr3": 5.8, "attr4": 320.0}}, {"type": "Feature", "id": 8, "geometry": {"type": "LineString", "coordinates": [[17, 34], [16.887493858644547, 34.13356039140768], [17.050774167231765, 34.269009846086995], [16.93790920922015, 34.40256384863036], [17.101712917701974, 34.53800624320647], [16.962186228592657, 34.67529960820302], [17.126527752866963, 34.81073481562158], [16.98654060783366, 34.948021207463654]]}, "properties": {"start-timestamp": "2021-04-02 12:00:00", "end-timestamp": "2021-04-02 19:00:00"}}]}

image

I think using getSpeed instead of getRoutageSpeed is correct, since creating isochrones we need the effective speed retrieved from polar (or an interpolated value for missing angles). getRoutageSpeed instead, gives a value I'm not quite sure what it is.

Please give it a try and let me know

enricofer commented 2 years ago

You are right. using getSpeed i'm getting much more reasonable results: https://gist.github.com/enricofer/e053c9ce45fb20d7cbafd595b17ab8de The getRoutageSpeed method as far as I understand calculate the best VGM https://it.wikipedia.org/wiki/Velocity_Made_Good for a gived wind relative angle, but iterating over ALL direction to find the most effective route is a concept that overcomes and contains the search for the the best VMG. Comparing the test results for getSpeed evaluation method with those for getRoutageSpeed evaluation method: https://gist.github.com/enricofer/411843c5d1cd9e73292d84ebf9880358, keeping apart the faulty 000 heading, appears that the methods outputs a longest route but allowing faster speed than the second method but the reaching time is the same, and therefore, since wind propellent is free, the result is the same, but increasing algorithm performances. For completeness, tests should be performed in a more real conditions but I'm expecting the same result so I would get rid of getRoutageSpeed method for best route evaluation. If you agree, I would submit a PR. ping @pcav

pcav commented 2 years ago

Fine to me, thanks for letting me know @enricofer

dakk commented 2 years ago

@enricofer agree for the PR, getSpeed is the correct approach for isochrones calculations: the VMG is intrinsic here, since the isochrones tree will return fastest route anyway. We can create another router class using VMG calculations without isochrones.

Waiting for PR.

dakk commented 2 years ago

Closed by PR #9