pretty plot #93

opendrift_leeway_webgui/simulation.py

@@ -11,6 +11,14 @@
 import os
 import uuid
 from datetime import datetime, timedelta
+import numpy as np
+from matplotlib import ticker
+import matplotlib.pyplot as plt
+from matplotlib.colors import ListedColormap
+from matplotlib.collections import LineCollection
+import cartopy.crs as ccrs
+from  cartopy.mpl import gridliner
+
 
 # pylint: disable=import-error
 from opendrift.models.leeway import Leeway
@@ -107,11 +115,177 @@ 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]):
+        if np.isnan(lon[i,-1]):
+            lon_max_idx=np.where(np.isnan(lon[i,...]))[0][0]-1
+        else:
+            lon_max_idx=-1
+        if np.isnan(lat[i,-1]):
+            lat_max_idx=np.where(np.isnan(lat[i,...]))[0][0]-1
+        else:
+            lat_max_idx=-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 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):
+    deg=np.floor(decimal)
+    min=(decimal-deg)*60
+    min=np.floor(min)
+    return deg+(min/60)
+
+def ceil_min(decimal):
+    deg=np.floor(decimal)
+    min=(decimal-deg)*60
+    min=np.ceil(min)
+    return deg+(min/60)
+
+def get_zebra_line(points,gcrs):
+    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()