diff --git a/PSILARTEMP.py b/PSILARTEMP.py index 57cbad1e1bf77d58f074e1dea5b8b98197f1935d..b3fd6c1ab72109fd2ffdcdd12f80b7b150c0bcbf 100755 --- a/PSILARTEMP.py +++ b/PSILARTEMP.py @@ -16,7 +16,8 @@ along with this program. If not, see <https://www.gnu.org/licenses/>. """ -import subprocess, os +import subprocess +import os from dotenv import load_dotenv from datetime import datetime import pytz @@ -30,6 +31,7 @@ load_dotenv() config = configparser.ConfigParser() config.read("PSILARTEMP.cfg") +year = os.getenv("YEAR") weatherdata_path = os.getenv("WEATHER_DATA_DIR") tmp_path = "tmp/" out_path = os.getenv("DATA_DIR") @@ -37,19 +39,32 @@ out_path = os.getenv("DATA_DIR") # TODO: Put timezone out in .env 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) +subprocess.run( + f"cdo -selname,TM -seldate,{year}-04-01T00:00:00,{year}-12-31T00:00:00 {weatherdata_path}met_1_0km_nordic-{year}.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) +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) +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) +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) +subprocess.run( + f"cdo -timcumsum {tmp_path}TM_minus_5_nozero.nc {tmp_path}DD_unmasked.nc", + shell=True, +) # Mask results using a CSV file with polygons @@ -57,7 +72,10 @@ subprocess.run(f"cdo -timcumsum {tmp_path}TM_minus_5_nozero.nc {tmp_path}DD_unma if os.getenv("MASK_FILE") is not None: mask_file = os.getenv("MASK_FILE") print(f"Applying mask file {mask_file} to result.nc") - subprocess.run(f'cdo -maskregion,{mask_file} {tmp_path}DD_unmasked.nc {tmp_path}DD.nc', shell=True) + subprocess.run( + f"cdo -maskregion,{mask_file} {tmp_path}DD_unmasked.nc {tmp_path}DD.nc", + shell=True, + ) else: os.rename(f"{tmp_path}DD_unmasked.nc", f"{tmp_path}DD.nc") @@ -68,66 +86,87 @@ else: # 4: 360 < DD <= 560 # Peak flight period # 5: 560 < DD # 1st generation flight period ended # (A>0)*(A<=260)*2 + (A>260)*(A<=360)*3 + (A>360)*(A<=560)*4 + (A>560)*5 -subprocess.run(f'cdo -aexpr,"WARNING_STATUS=(TM>0)*(TM<=260)*2 + (TM>260)*(TM<=360)*3 + (TM>360)*(TM<=560)*4 + (TM>560)*5" {tmp_path}DD.nc {tmp_path}result.nc', shell=True) +subprocess.run( + f'cdo -aexpr,"WARNING_STATUS=(TM>0)*(TM<=260)*2 + (TM>260)*(TM<=360)*3 + (TM>360)*(TM<=560)*4 + (TM>560)*5" {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') +result = nc.Dataset(f"{tmp_path}result.nc", "r") timesteps = result.variables["time"][:] -#print(timesteps) +# 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}") + # 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"rm {tmp_path}result_*{file_date}_lcc.tif", shell=True) + subprocess.run(f"rm {out_path}result_*{file_date}.tif", shell=True) with open("/dev/null", "w") as devnull: - 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, stdout=devnull) - 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, stdout=devnull) + 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, + stdout=devnull, + ) + 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, + stdout=devnull, + ) # 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, stdout=devnull) - subprocess.run(f'gdalwarp -t_srs EPSG:4326 {tmp_path}result_{file_date}_lcc.tif {out_path}result_{file_date}.tif', shell=True, stdout=devnull) - #""" + 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, + stdout=devnull, + ) + subprocess.run( + f"gdalwarp -t_srs EPSG:4326 {tmp_path}result_{file_date}_lcc.tif {out_path}result_{file_date}.tif", + shell=True, + stdout=devnull, + ) + # """ 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(timestep_dates) -#print(len(timestep_dates)) +# print(len(timestep_dates)) # Generate mapfile # Building data sets for language specific legends languages = [] -language_codes = config["i18n"]["languages"].split(","); +language_codes = config["i18n"]["languages"].split(",") for language_code in language_codes: - language = {"language_code": language_code} - if ("i18n.%s" % language_code) in config: + language = {"language_code": language_code} + if ("i18n.%s" % language_code) in config: for keyword in config["i18n.%s" % language_code]: language[keyword] = config["i18n.%s" % language_code][keyword] languages.append(language) # The paths should be set in a .env file -env = Environment(loader=FileSystemLoader('.')) +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"), - "languages": languages, - "language_codes": language_codes -}) +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"), + "languages": languages, + "language_codes": language_codes, + } +) 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 +with open(f"{mapfile_outdir}/PSILARTEMP.map", "w") as f: + f.write(output) diff --git a/env-sample b/env-sample index f2617583b5b7abd5b8b5f2b82d60f8cf3b42d417..03888741aadefd402ce984960ddf29f7f479cc30 100644 --- a/env-sample +++ b/env-sample @@ -1,5 +1,7 @@ # Use this as an example to create your own .env file +# Year to run for +YEAR=2024 # Where your application resides HOME_DIR=/home/foo/ # Path to the weather data