Initial commit

.gitignore

+*.nc
+*.tif
+*.tif.aux.xml
+.env
+.venv
+mapfile/PSILARTEMP.map
\ No newline at end of file

PSILARTEMP.py

+#!/usr/bin/python3
+"""
+    Copyright (C) 2023 NIBIO <https://www.nibio.no/>. 
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Affero General Public License as
+    published by the Free Software Foundation, either version 3 of the
+    License, or (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Affero General Public License for more details.
+
+    You should have received a copy of the GNU Affero General Public License
+    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+"""
+
+import subprocess, os
+from dotenv import load_dotenv
+from datetime import datetime
+import pytz
+import netCDF4 as nc
+from jinja2 import Environment, FileSystemLoader
+
+load_dotenv()
+
+weatherdata_path="in/"
+tmp_path="tmp/"
+out_path="out/"
+
+# TODO: Make this truly independent of time zones/Norwegian conditions?
+utc_offset = "+02:00"
+local_timezone = pytz.timezone("Europe/Oslo")
+
+"""
+
+# Calculate cumulated degree days above 5 degrees after 1st of April
+# Remove all values before April 1st
+subprocess.run(f"cdo -selname,TM -seldate,2023-04-01T00:00:00,2023-12-31T00:00:00 {weatherdata_path}met_1_0km_nordic-2023.nc {tmp_path}TM_from_april.nc", shell=True)
+# Subtracting 5 deg C from all cells
+subprocess.run(f"cdo -subc,5 {tmp_path}TM_from_april.nc {tmp_path}TM_minus_5.nc", shell=True)
+# Create an .nc file with all cells containing 1 if cell in tgminus5.nc is >=0, 0 otherwise
+subprocess.run(f"cdo -gtc,0 {tmp_path}TM_minus_5.nc {tmp_path}TM_minus_5_gtc.nc", shell=True)
+
+# Multiplying tgminus5.nc with tgminus5gtc.nc, so that all sub zero values are set to 0
+subprocess.run(f"cdo -mul {tmp_path}TM_minus_5.nc {tmp_path}TM_minus_5_gtc.nc {tmp_path}TM_minus_5_nozero.nc", shell=True)
+# Accumulate day degrees with the corrected values
+subprocess.run(f"cdo -timcumsum {tmp_path}TM_minus_5_nozero.nc {tmp_path}DD_unmasked.nc", shell=True)
+
+
+
+# Mask the output with Norway land borders
+subprocess.run(f"cdo -maskregion,Norge_landomrader.csv {tmp_path}DD_unmasked.nc {tmp_path}DD.nc", shell=True)
+
+# Add the DD threshold classification => warning status
+# (A>0)*(A<=260)*2 + (A>260)*(A<=360)*3 + (A>360)*(A<=560)*4 + (A>560)*0
+subprocess.run(f'cdo -aexpr,"WARNING_STATUS=(TM>0)*(TM<=260)*2 + (TM>260)*(TM<=360)*3 + (TM>360)*(TM<=560)*4 + (TM>560)*0" {tmp_path}DD.nc {tmp_path}result.nc', shell=True)
+
+"""
+
+# Split the combined file into daily .nc files again, with YYYY-MM-DD in the filename. Convert to corresponding GeoTIFF files
+# Variables that needs discrete classification, must be integers in order for mapserver to work properly (Don't ask!)
+# Since we need WARNING_STATUS to be discretely classified, we need to create a separate GeoTIFF file for it
+result = nc.Dataset(f'{tmp_path}result.nc', 'r')
+timesteps = result.variables["time"][:]
+#print(timesteps)
+timestep_index = 1
+timestep_dates = []
+for timestep in timesteps:
+    timestep_date = datetime.fromtimestamp(timestep).astimezone(local_timezone)
+    file_date = timestep_date.astimezone(local_timezone).strftime("%Y-%m-%d")
+    print(f"{timestep_index}: {file_date}")
+    timestep_dates.append(file_date)
+    # We must delete the GeoTIFF file before merging
+    """subprocess.run(f'rm {tmp_path}result_*{file_date}_lcc.tif', shell=True)
+    subprocess.run(f'rm {out_path}result_*{file_date}.tif', shell=True)
+    subprocess.run(f'gdal_translate -ot Int16 -of GTiff  -b {timestep_index} NETCDF:"{tmp_path}result.nc":WARNING_STATUS {tmp_path}result_WARNING_STATUS_{file_date}_lcc.tif', shell=True)
+    subprocess.run(f'gdal_translate -ot Float32 -of GTiff -b {timestep_index} NETCDF:"{tmp_path}result.nc":TM {tmp_path}result_{file_date}_lcc.tif', shell=True)
+    # Need to reproject the files, to ensure we have the projection given in the generted mapfile. We always use EPSG:4326 for this
+    subprocess.run(f'gdalwarp -t_srs EPSG:4326 {tmp_path}result_WARNING_STATUS_{file_date}_lcc.tif {out_path}result_WARNING_STATUS_{file_date}.tif', shell=True)
+    subprocess.run(f'gdalwarp -t_srs EPSG:4326 {tmp_path}result_{file_date}_lcc.tif {out_path}result_{file_date}.tif', shell=True)
+    """
+    timestep_index = timestep_index + 1
+
+if len(timestep_dates) != len(set(timestep_dates)):
+    print("ERROR: Something is wrong with your timesteps")
+    #print(timestep_dates)
+
+print(len(timestep_dates))
+
+# Generate mapfile
+# The paths should be set in a .env file
+env = Environment(loader=FileSystemLoader('.'))
+template = env.get_template("mapfile/template.j2")
+output = template.render({
+     "timestep_dates": timestep_dates,
+     "mapserver_data_dir": os.getenv("MAPSERVER_DATA_DIR"),
+     "mapserver_mapfile_dir": os.getenv("MAPSERVER_MAPFILE_DIR"),
+     "mapserver_log_file": os.getenv("MAPSERVER_LOG_FILE"),
+     "mapserver_image_path": os.getenv("MAPSERVER_IMAGE_PATH"),
+     "mapserver_extent": os.getenv("MAPSERVER_EXTENT")
+})
+mapfile_outdir = os.getenv("MAPFILE_DIR")
+with open(f"{mapfile_outdir}/PSILARTEMP.map", 'w') as f:
+        f.write(output)
\ No newline at end of file

in/.gitignore

+*
+*/
+!.gitignore
\ No newline at end of file

mapfile/query_template.xml

+<!--mapserver template-->
+<?xml version="1.0" encoding="UTF-8"?>
+<vipsResult>
+    <modelName value="Gulrotflue svermetidspunktmodell"/>
+    <modelId value="PSILARTEMP"/>
+    <warningStatus value="[value_0]"/>
+</vipsResult>
\ No newline at end of file

mapfile/query_template_DD.xml

+<!--mapserver template-->
+<?xml version="1.0" encoding="UTF-8"?>
+<vipsResult>
+    <modelName value="Gulrotflue svermetidspunktmodell"/>
+    <modelId value="PSILARTEMP"/>
+    <parameter name="DD" value="[value_0]"/>
+</vipsResult>
\ No newline at end of file

mapfile/template.j2

+MAP
+
+# This mapfile is generated using Jinja2 templates
+# 
+# Copyright (C) 2023 NIBIO <https://www.nibio.no/>
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+NAME "VIPS.PSILARTEMP"
+
+# WGS84
+EXTENT {{mapserver_extent}}
+UNITS DD
+
+CONFIG "MS_ERRORFILE" "{{mapserver_log_file}}" 
+
+PROJECTION 
+    "init=epsg:4326"
+END
+
+
+WEB
+    IMAGEPATH "{{mapserver_image_path}}"
+    IMAGEURL "/tmp/"
+    # List of standard metadata: https://mapserver.org/ogc/wms_server.html#web-object-metadata
+    # i18n support: https://mapserver.org/ogc/inspire.html#inspire-multi-language-support  
+    METADATA
+            "wms_keywordlist" "VIPS model Carrot Rust Fly (PSILARTEMP)"
+            "wms_abstract"  "
+                <p>The warning system model &laquo;Carrot rust fly temperature&raquo; is based on a Finnish temperature-based model (Markkula <em>et al</em>, 1998; Tiilikkala &amp; Ojanen, 1999; Markkula <em>et al</em>, 2000). The model determines the start of the flight period for the 1st and 2nd generation of carrot rust fly based on accumuleted degree-days (day-degrees) over a base temperature of 5,0 &deg;C. VIPS uses the model for the 1st generation only.</p>
+                <p>Standard air temperature (temperature measured 2 m above ground) is used in the model. Degree-days are defined for this model as the sum of the difference between a base temperature of 5,0 &deg;C and the mean temperature for all days with a temperature &gt;5,0 &deg;C, in other words (daily mean temperature &ndash; 5,0 &deg;C) from 1 March (beginning when the ground has thawed).</p>
+                <h3>Interpretation of the warning</h3>
+                <p>Green rectangles indicate that the flight period has not yet begun. The accumulated day-degrees are &lt; 260 degree-days (day-degrees).</p>
+                <p>Yellow rectangles indicate that the flight period is beginning and that flies can be coming into the field. The accumulated day-degrees are &ge; 260 degree-days (day-degrees).</p>
+                <p>Red rectangles indicate peak flight period. The accumulated day-degrees are &ge; 360 degree-days (day-degrees).</p>
+                <p>Grey rectangles indicate that the flight period of the 1st generation is over. The accumulated day-degrees are &ge; 560 degree-days (day-degrees).</p>
+                <p>Be aware that in areas with field covers (plastic, single or double non-woven covers, etc.) with early crops the preceding season (either on the current field or neighboring fields), the flight period can start earlier due to higher soil temperature under the covers.</p>
+                <p>The graphic presentation shows the progress of the day-degree calculation that is the basis for the green, yellow and red warnings (the line: &laquo;Accumulated day-degrees with the base temperature of 5 degrees C&raquo;). The graph show straight horizontal lines for the threshold values. &laquo;Day-degree threshold for low likelihood of attack&raquo; corresponds to when the warning changes from green to yellow (260 day-degrees). &laquo;Day-degree threshold for high likelihood of attack &raquo; corresponds to when the warning changes from yellow to red (360 day-degrees). &laquo;Day-degree threshold for end of flight period&raquo; corresponds to when the warning changes from red to grey (560 day-degrees) and the flight period for the 1st generation is considered to be over. When the line &laquo; Accumulated day-degrees with the base temperature of 5 degrees C&raquo; intersects one of the lines for the day-degree threshold, the warning will advance to the next level and the color of the warning rectangle will change. The graph also shows the daily mean temperature for the relevant weather station. The graph is dynamic and the user can choose which parameters that are shown by clicking on the explanation below the graph.</p>
+                <h3>Warning season &ndash; start and end of the warning</h3>
+                <p>Starting time: Shown in VIPS from 1 April (starting date for accumulation of degree-days beginning from when the ground has thawed).<br /> Ending time: When the model has reached the requirement for the end of the flight period for the 1st generation (560 day-degrees).</p>
+                <h3>Testing and validation of the model</h3>
+                <h4>National</h4>
+                <p>The model has been tested for Norwegian conditions in the period of 2011-2015. Validation has shown that the model was accurate for most locations for the 1st generation. The model was, however, inaccurate for the 2nd generation for many locations and was totally incorrect for the district of J&aelig;ren. Based on the validation results, we have therefore chosen to remove warning of the 2nd generation of carrot rust fly from this model.</p>
+                <h4>International</h4>
+                <p>The model is in use in Finland. The extent to which the model has been validated in Finland is uncertain, but a validation was done in the summer of 1997. This validation showed that the model was quite accurate, but that the threshold temperature sm should be lowered to 255 day-degrees for the start of the flight period (as opposed to 260 degree-days in the original model) and 355 degree-days forthe peak flight period (as opposed to 360 degree-days in the original model) (Tiilikkala &amp; Ojanen, 1999).</p>
+                <h3>Literature</h3>
+                <p>Markkula, I., H. Ojanen and K. Tiilikkala. 1998. Forecasting and monitoring of the carrot fly (<em>Psila rosae</em>) in Finland. The 1998 Brighton Conference: Pests and Diseases: Conference Proceedings, Volume 2, pages 657-662.</p>
+                <p>Tiilikkala, K. and H. Ojanen. 1999. Use of a geographical information system (GIS) for forecasting the activities of carrot fly and cabbage root fly. IOBC/WPRS Bulletin 22 (5): 15-24.</p>
+                <p>Markkula, I., A. Hannukkala, A. Lehtinen, I. Mattila, H. Ojanen, S. Raisjio, P. Reijonen and K. Tiilikkala. 2000. Pest warnings and forecasts sent as SMS messages from models into \"farmer's pocket\". In publication: Pests &amp; Diseases 2000. Proceedings of an international conference held at the Brighton Hilton Metropole Hotel, UK, 13-16 November 2000. BCPC Conference Proceedings pages 285-290.</p>
+                <p>&nbsp;</p>
+                <p><strong>Contacts</strong></p>
+                <p>Annette Folkedal Schj&oslash;ll <a href=\"mailto:annette.folkedal.schjoll@nibio.no\">annette.folkedal.schjoll@nibio.no</a><br />Tor J. Johansen <a href=\"mailto:tor.johansen@nibio.no\">tor.johansen@nibio.no</a></p>
+            "
+        "wms_enable_request"  "*"
+        "wms_title"     "Carrot rust fly (Psila rosae) temperature model"
+        "wms_getfeatureinfo_formatlist" "text/plain,text/html,text/xml"
+        "wms_accessconstraints" "none"
+        "wms_addresstype" ""
+        "wms_address" "Høgskoleveien 7"
+        "wms_city" "Ås"
+        "wms_stateorprovince" "Ås"
+        "wms_postcode" "1430"
+        "wms_country" "Norway"
+        "wms_contactelectronicmailaddress" "vips@nibio.no"
+        "wms_feature_info_mime_type" "text/html"
+        "wms_contactperson" "Berit Nordskog"
+        "wms_contactposition" "owner"
+        "wms_contactorganization" "Norsk institutt for bioøkonomi (NIBIO)"
+    END
+END #web
+
+############################# Start of legend #########################################################
+# The legend configuration uses the CLASSITEM information in each layer to generate colour icons and labels
+LEGEND
+    STATUS ON
+    KEYSIZE 16 12
+    LABEL
+        TYPE TRUETYPE
+        SIZE 10
+        COLOR 0 0 0
+        OFFSET 0 -3
+    END
+END
+
+{% for timestep_date in timestep_dates %}
+LAYER
+    NAME "PSILARTEMP.WARNING_STATUS.{{ timestep_date }}"
+    DATA "{{mapserver_data_dir}}result_WARNING_STATUS_{{ timestep_date }}.tif"
+    TEMPLATE "{{mapserver_mapfile_dir}}query_template.xml" TOLERANCE 1 TOLERANCEUNITS PIXELS
+    TYPE RASTER
+    PROCESSING "BANDS=1" # WARNING_STATUS band on top (others invisible, but band values are available in the query template)
+    PROCESSING "NODATA=-1"
+   
+    
+    STATUS ON
+    METADATA
+      "wms_title"     "Carrot rust fly (Psila rosae) temperature model {{ timestep_date }}"  
+    END  
+      CLASSITEM "[pixel]"
+  
+    # class using simple string comparison, equivalent to ([pixel] = 0)
+  
+    CLASS
+      NAME "Model not running"
+      EXPRESSION ([pixel] >= 0 AND [pixel] < 2) 
+      STYLE
+          COLOR 112 112 112
+      END
+    END
+    CLASS
+      NAME "No infection risk"
+      EXPRESSION ([pixel] >= 2 AND [pixel] < 3) 
+      STYLE
+          COLOR 0 180 87
+      END
+    END
+    CLASS
+      NAME "Possible infection risk"
+      EXPRESSION ([pixel] >= 3 AND [pixel] < 4) 
+      STYLE
+          COLOR 255 204 0
+      END
+    END
+    CLASS
+      NAME "High infection risk"
+      EXPRESSION ([pixel] >= 4) 
+      STYLE
+          COLOR  255 0 0
+      END
+    END
+END # Layer
+  
+LAYER
+     NAME "PSILARTEMP.DD.{{ timestep_date }}"
+     DATA "{{mapserver_data_dir}}result_{{ timestep_date }}.tif"
+     TEMPLATE "{{mapserver_mapfile_dir}}query_template_DD.xml" TOLERANCE 1 TOLERANCEUNITS PIXELS
+     TYPE RASTER
+     #PROCESSING "BANDS=1" # DD band on top (others invisible, but band values are available in the query template)
+     #PROCESSING "SCALE=AUTO"
+     #PROCESSING "NODATA=-1"
+   
+   
+     STATUS ON
+     METADATA
+       "wms_title"     "Carrot rust fly (Psila rosae) temperature model Day Degrees {{ timestep_date }}"
+     END
+    CLASSITEM "[pixel]"
+    CLASS
+      NAME "Day degrees"
+      #EXPRESSION ([pixel] >= 0 AND [pixel] <= 72)
+      STYLE
+        DATARANGE 0 1000
+        COLORRANGE 0 0 255 255 0 0
+      END
+    END
+   
+END # Layer
+
+{% endfor %}
+  
+END #map
+  
+  
\ No newline at end of file

out/.gitignore

+*
+*/
+!.gitignore
\ No newline at end of file

tmp/.gitignore

+*
+*/
+!.gitignore
\ No newline at end of file