Plot simulation results on nautical map, fix #93

.github/workflows/linting.yml

@@ -23,7 +23,7 @@ jobs:
       - uses: actions/setup-python@v4
       - uses: jamescurtin/isort-action@master
         with:
-          configuration: --check-only
+          configuration: --check-only -n
   pylint:
     runs-on: ubuntu-latest
     steps:

.github/workflows/tests.yml

@@ -19,5 +19,16 @@ jobs:
         uses: supercharge/[email protected]
       - name: Install dependencies
         run: pip install -e .[dev]
+      - name: Login to Docker Hub
+        uses: docker/login-action@v3
+        with:
+          username: ${{ secrets.DOCKERHUB_USERNAME }}
+          password: ${{ secrets.DOCKERHUB_TOKEN }}
+      - name: Build custom container with OWSLib
+        uses: docker/build-push-action@v5
+        with:
+          file: opendrift/Dockerfile
+          push: false
+          tags: opendrift-leeway-custom:latest
       - name: Run tests
         run: pytest

README.md

@@ -57,6 +57,11 @@ Authentication can be provided in two ways:
    python3 manage.py migrate
    python3 manage.py createsuperuser
    ```
+6. Build the enhanced Docker container:
+   ```bash
+   cd opendrift
+   docker build -t opendrift-leeway-custom .
+   ```
 
 # Development Server
 

opendrift/Dockerfile

@@ -0,0 +1,2 @@
+FROM opendrift/opendrift
+RUN pip3 install OWSLib

opendrift_leeway_webgui/api/urls.py

@@ -1,6 +1,7 @@
 """
 URL patterns for the Opendrift Leeway Webgui API
 """
+
 from django.urls import include, path
 
 #: The namespace for this URL config (see :attr:`django.urls.ResolverMatch.app_name`)

opendrift_leeway_webgui/api/v1/urls.py

@@ -1,6 +1,7 @@
 """
 URL patterns for the first version of the Opendrift Leeway Webgui API
 """
+
 from django.urls import include, path
 from drf_spectacular.views import SpectacularAPIView, SpectacularSwaggerView
 from rest_framework.routers import DefaultRouter

opendrift_leeway_webgui/core/settings.py

@@ -7,6 +7,7 @@
 For the full list of settings and their values, see
 https://docs.djangoproject.com/en/3.2/ref/settings/
 """
+
 import os
 from pathlib import Path
 

opendrift_leeway_webgui/core/urls.py

@@ -13,6 +13,7 @@
     1. Import the include() function: from django.urls import include, path
     2. Add a URL to urlpatterns:  path('blog/', include('blog.urls'))
 """
+
 from django.contrib import admin
 from django.urls import include, path
 

opendrift_leeway_webgui/leeway/celery.py

@@ -1,6 +1,7 @@
 """
 Configure Celery workers
 """
+
 import configparser
 import email
 import os

opendrift_leeway_webgui/leeway/forms.py

@@ -1,6 +1,7 @@
 """
 Forms for the web GUI
 """
+
 from django.forms import CharField, ModelForm, TextInput
 
 from .models import LeewaySimulation

opendrift_leeway_webgui/leeway/tasks.py

@@ -31,7 +31,7 @@ def run_leeway_simulation(request_id):
         f"{settings.SIMULATION_ROOT}:/code/leeway",
         "--volume",
         f"{settings.SIMULATION_SCRIPT_PATH}:/code/leeway/simulation.py",
-        "opendrift/opendrift",
+        "opendrift-leeway-custom:latest",
         "python3",
         "leeway/simulation.py",
         "--longitude",

opendrift_leeway_webgui/leeway/urls.py

@@ -4,6 +4,7 @@
 
 For more information on this file, see :doc:`django:topics/http/urls`.
 """
+
 from django.conf import settings
 from django.conf.urls.static import static
 from django.urls import include, path

opendrift_leeway_webgui/leeway/utils.py

@@ -1,6 +1,7 @@
 """
 Utilities
 """
+
 from django.apps import apps
 from django.contrib.auth import get_user_model
 from django.core.mail import EmailMessage, send_mail

opendrift_leeway_webgui/simulation.py

@@ -1,3 +1,4 @@
+# pylint: disable=too-many-function-args
 """
 Wrapper/Helper script for Leeway simulations with OpenDrift.
 
@@ -13,12 +14,24 @@
 from datetime import datetime, timedelta
 
 # pylint: disable=import-error
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+import numpy as np
+
+# pylint: disable=import-error
+from cartopy.mpl import gridliner
+from matplotlib import ticker
+from matplotlib.collections import LineCollection
+from matplotlib.colors import ListedColormap
+
+# pylint: disable=import-error, disable=no-name-in-module
 from opendrift.models.leeway import Leeway
 from opendrift.readers import reader_global_landmask
 
 INPUTDIR = "/code/leeway/input"
 
 
+# pylint: disable=too-many-locals, disable=too-many-statements
 def main():
     """Run opendrift leeway simulation"""
     parser = argparse.ArgumentParser(description="Simulate drift of object")
@@ -107,11 +120,229 @@ def main():
     simulation.run(
         duration=timedelta(hours=args.duration), time_step=600, outfile=f"{outfile}.nc"
     )
-    simulation.plot(
-        fast=True, legend=True, filename=f"{outfile}.png", linecolor="age_seconds"
+
+    # Plotting results
+    lon, lat = np.array(simulation.get_lonlats())
+    lon[lon == 0] = np.nan
+    lat[lat == 0] = np.nan
+
+    crs = ccrs.Mercator()  # Mercator projection to have angle true projection
+    gcrs = ccrs.PlateCarree(globe=crs.globe)  # PlateCarree for straight lines
+
+    fig = plt.figure(figsize=(8, 8))  # figsize set low to get small files
+    ax = plt.axes(projection=crs)
+
+    # base map layer
+    ax.add_wms(wms="https://sgx.geodatenzentrum.de/wms_topplus_open", layers=["web"])
+    # quote source: Kartendarstellung: © Bundesamt für Kartographie und Geodäsie
+    # (2021), Datenquellen:
+    # https://gdz.bkg.bund.de/index.php/default/wms-topplusopen-wms-topplus-open.html
+
+    stranded = None
+    active = None
+    for i in range(lon.shape[0]):
+        lon_max_idx = (
+            np.where(np.isnan(lon[i, ...]))[0][0] - 1 if np.isnan(lon[i, -1]) else -1
+        )
+        lat_max_idx = (
+            np.where(np.isnan(lat[i, ...]))[0][0] - 1 if np.isnan(lat[i, -1]) else -1
+        )
+
+        points = np.array([lon[i, ...], lat[i, ...]]).T.reshape(-1, 1, 2)
+        segments = np.concatenate([points[:-1], points[1:]], axis=1)
+
+        lc = LineCollection(
+            segments, cmap="jet", norm=plt.Normalize(0, args.duration), transform=gcrs
+        )
+        lc.set_array(np.linspace(0, args.duration, len(segments)))
+        lc.set_linewidth(2)
+        line = ax.add_collection(lc)
+
+        initial = ax.scatter(
+            lon[i, 0],
+            lat[i, 0],
+            transform=gcrs,
+            color="green",
+            zorder=100 + i,
+            s=15,
+            edgecolor="black",
+            linewidth=0.5,
+        )
+        if lon_max_idx > 0:
+            stranded = ax.scatter(
+                lon[i, lon_max_idx],
+                lat[i, lat_max_idx],
+                transform=gcrs,
+                color="red",
+                zorder=100 + i + 1,
+                s=15,
+                edgecolor="black",
+                linewidth=0.5,
+            )
+        else:
+            active = ax.scatter(
+                lon[i, lon_max_idx],
+                lat[i, lat_max_idx],
+                transform=gcrs,
+                color="blue",
+                zorder=100 + i + 1,
+                s=15,
+                edgecolor="black",
+                linewidth=0.5,
+            )
+
+    # print legend for points
+    initial.set_label("Initial")
+    if stranded is not None:
+        stranded.set_label("Stranded")
+    if active is not None:
+        active.set_label("Active")
+
+    ax.legend(loc="center right", bbox_to_anchor=(0, 0.5))
+
+    # add colorbar with hours
+    cbar = fig.colorbar(line, location="bottom")
+    cbar.set_label("Hours")
+
+    # set map extent
+    x_max = np.nanmax(lon)
+    x_min = np.nanmin(lon)
+    y_max = np.nanmax(lat)
+    y_min = np.nanmin(lat)
+
+    dis_x = (x_max - x_min) * 0.02
+    dis_y = (y_max - y_min) * 0.02
+    extent = [
+        floor_min(x_min - dis_x),
+        ceil_min(x_max + dis_x),
+        floor_min(y_min - dis_y),
+        ceil_min(y_max + dis_y),
+    ]
+    ax.set_extent(extent)
+
+    # prepare gridlines on good position
+    x_step = 1 / 60  # step size for grid lines a line every minute
+    x_step_div = 10  # num of zebra stripes between grid lines
+    y_step = 1 / 60
+    y_step_div = 10
+
+    x_ticks = np.arange(extent[0], extent[1], x_step)
+    y_ticks = np.arange(extent[2], extent[3], y_step)
+
+    # check if too many grid lines:
+    if len(x_ticks) > 100:  # a line every 5 minutes
+        x_step = 1 / 12
+        x_step_div = 5
+        x_ticks = np.arange(extent[0], extent[1], x_step)
+    if len(y_ticks) > 100:
+        y_step = 1 / 12
+        y_step_div = 5
+        y_ticks = np.arange(extent[2], extent[3], y_step)
+    if len(x_ticks) > 100:  # a line every 10 minutes
+        x_step = 1 / 6
+        x_step_div = 10
+        x_ticks = np.arange(extent[0], extent[1], x_step)
+    if len(y_ticks) > 100:
+        y_step = 1 / 6
+        y_step_div = 10
+        y_ticks = np.arange(extent[2], extent[3], y_step)
+
+    xloc = ticker.FixedLocator(x_ticks)
+    yloc = ticker.FixedLocator(y_ticks)
+
+    longitude_formatter = gridliner.LongitudeFormatter(dms=True, auto_hide=False)
+    latitude_formatter = gridliner.LatitudeFormatter(dms=True, auto_hide=False)
+
+    gl = ax.gridlines(
+        draw_labels=True, linewidth=1, color="gray", alpha=0.5, rotate_labels=True
+    )
+    gl.ylocator = yloc
+    gl.xlocator = xloc
+    gl.xlabels_top = False
+    gl.ylabels_left = False
+    gl.xformatter = longitude_formatter
+    gl.yformatter = latitude_formatter
+
+    # add zebra frame
+    zebra_x = np.arange(extent[0], extent[1] + x_step / x_step_div, x_step / x_step_div)
+    zebra_y = np.arange(extent[2], extent[3] + y_step / y_step_div, y_step / y_step_div)
+
+    if len(zebra_x) > 200 or len(zebra_y) > 200:
+        x_step_div = 2
+        zebra_x = np.arange(
+            extent[0], extent[1] + x_step / x_step_div, x_step / x_step_div
+        )
+        y_step_div = 2
+        zebra_y = np.arange(
+            extent[2], extent[3] + y_step / y_step_div, y_step / y_step_div
+        )
+
+    points = np.array([zebra_x, np.zeros_like(zebra_x) + extent[2]]).T.reshape(-1, 1, 2)
+    lc = get_zebra_line(points, gcrs)
+    line = ax.add_collection(lc)
+    points = np.array([zebra_x, np.zeros_like(zebra_x) + extent[3]]).T.reshape(-1, 1, 2)
+    lc = get_zebra_line(points, gcrs)
+    line = ax.add_collection(lc)
+    points = np.array([np.zeros_like(zebra_y) + extent[0], zebra_y]).T.reshape(-1, 1, 2)
+    lc = get_zebra_line(points, gcrs)
+    line = ax.add_collection(lc)
+    points = np.array([np.zeros_like(zebra_y) + extent[1], zebra_y]).T.reshape(-1, 1, 2)
+    lc = get_zebra_line(points, gcrs)
+    line = ax.add_collection(lc)
+
+    start = datetime.strptime(args.start_time, "%Y-%m-%d %H:%M")
+    end = start + timedelta(hours=args.duration)
+
+    plt.title(
+        f"Leeway Simulation Object Type: {simulation.get_config('seed:object_type')}\n"
+        f" From {start.strftime('%d-%m-%Y %H:%M')} to {end.strftime('%d-%m-%Y %H:%M')} UTC"
+    )
+    fig.text(
+        0,
+        0,
+        "Kartendarstellung: © Bundesamt für Kartographie und Geodäsie (2021),\nDatenquellen:"
+        " https://gdz.bkg.bund.de/index.php/default/wms-topplusopen-wms-topplus-open.html",
+        fontsize=8,
     )
+
+    fig.savefig(f"{outfile}.png")
+
     print(f"Success: {outfile}.png written.")
 
 
+def floor_min(decimal):
+    """
+    floor funtion for arc minutes
+    """
+    deg = np.floor(decimal)
+    minutes = (decimal - deg) * 60
+    minutes = np.floor(minutes)
+    return deg + (minutes / 60)
+
+
+def ceil_min(decimal):
+    """
+    ceil fnction for arc minutes
+    """
+    deg = np.floor(decimal)
+    minutes = (decimal - deg) * 60
+    minutes = np.ceil(minutes)
+    return deg + (minutes / 60)
+
+
+def get_zebra_line(points, gcrs):
+    """
+    Define a blak/white dashed line for a map frame
+    """
+    segments = np.concatenate([points[:-1], points[1:]], axis=1)
+    cmap = ListedColormap(["k", "w"])
+    lc = LineCollection(segments, cmap=cmap, norm=plt.Normalize(0, 1), transform=gcrs)
+    array = np.zeros(points.shape[0])
+    array[::2] = 1
+    lc.set_array(array)
+    lc.set_linewidth(4)
+    return lc
+
+
 if __name__ == "__main__":
     main()

pyproject.toml

@@ -36,6 +36,8 @@ dependencies = [
 dev = [
     "bumpver",
     "pre-commit",
+    "numpy",
+    "matplotlib",
     "pylint",
     "pylint-django",
     "pytest",

tests/conftest.py

@@ -1,6 +1,7 @@
 """
 This module contains shared fixtures for pytest
 """
+
 import pytest
 
 from .leeway.test_simulation import _teardown_test_simulation