diff --git a/PSILARTEMP.py b/PSILARTEMP.py index 315e0465986a0c713f8ebcb058228b995a196f4f..7bb0a2215b6ff46a0cd815e0da5a78944daeb5b2 100755 --- a/PSILARTEMP.py +++ b/PSILARTEMP.py @@ -48,10 +48,14 @@ subprocess.run(f"cdo -mul {tmp_path}TM_minus_5.nc {tmp_path}TM_minus_5_gtc.nc {t 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 -# TODO make optional -subprocess.run(f"cdo -maskregion,Norge_landomrader.csv {tmp_path}DD_unmasked.nc {tmp_path}DD.nc", shell=True) +# Mask results using a CSV file with polygons +# Env variable MASK_FILE must be set +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) +else: + os.rename(f"{tmp_path}DD_unmasked.nc", f"{tmp_path}DD.nc") # 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 @@ -68,23 +72,24 @@ 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'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) + 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) + # 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) 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)) +#print(len(timestep_dates)) # Generate mapfile # The paths should be set in a .env file