Commit 2f3de93d authored by Simon Gascoin's avatar Simon Gascoin
Browse files

Merge branch...

Merge branch '25-add-tools-to-compute-psl-and-nobs-qc-layer-of-psl-from-theia-products' into 'develop'

Resolve "Add tools to compute PSL and Nobs (QC layer of PSL) from Theia products"

Closes #25

See merge request !4
parents 03eaf057 eb20c702
#!/bin/bash
#PBS -N TheiaNeige
#PBS -J 1-191:1
#PBS -l select=1:ncpus=1:mem=20000mb
#PBS -l walltime=00:59:00
#PBS -M gascoins@cesbio.cnes.fr
#PBS -m e
# Load lis environnment
module load lis/develop
# Load all the available product names from the tile directory
pin=/work/OT/siaa/Theia/Neige/SNOW_ANNUAL_MAP_LIS_1.5/S2_with_L8_Densification/
inputFiles=($(find $pin -name multitemp_cloud_mask.vrt))
# use the PBS_ARRAY_INDEX variable to distribute jobs in parallel (bash indexing is zero-based)
i="${inputFiles[${PBS_ARRAY_INDEX} - 1]}"
# run script
python "${PBS_O_WORKDIR}"/compute_NOBS.py ${i}
#!/bin/bash
#PBS -N TheiaNeige
#PBS -J 1-191:1
#PBS -l select=1:ncpus=1:mem=20000mb
#PBS -l walltime=00:59:00
#PBS -M gascoins@cesbio.cnes.fr
#PBS -m e
# Load lis environnment
module load lis/develop
# Load all the available product names from the tile directory
pin=/work/OT/siaa/Theia/Neige/SNOW_ANNUAL_MAP_LIS_1.5/S2_with_L8_Densification/
inputFiles=($(find $pin -maxdepth 3 -name multitemp_cloud_mask.vrt))
# use the PBS_ARRAY_INDEX variable to distribute jobs in parallel (bash indexing is zero-based)
i="${inputFiles[${PBS_ARRAY_INDEX} - 1]}"
# run script
python "${PBS_O_WORKDIR}"/compute_PSL.py ${i}
#!/usr/bin/env python
# This script computes NOBS.tif, the number of clear observations to compute the SCD, SMOD and SOD syntheses
# In HAL the dependencies are loaded with module load lis/develop
# Author: Simon Gascoin
import rasterio
import numpy as np
import os,sys
# input file is the cloud mask vrt
# Example: /work/OT/siaa/Theia/Neige/SNOW_ANNUAL_MAP_LIS_1.5/S2_with_L8_Densification//T32TLP_20180901_20190831/tmpdir/multitemp_cloud_mask.vrt
f=sys.argv[1]
src=rasterio.open(f, 'r')
n=src.meta["count"]
W = src.read(range(1,n))
S=n-np.sum(W,axis=0)
outdir=os.path.split(os.path.split(f)[0])[0]
outfile=os.path.split(outdir)[1]
# export NOBS.tif in the parent folder of the input file
with rasterio.Env():
profile = src.profile
profile.update(
dtype=rasterio.uint16,
driver='GTiff',
count=1)
with rasterio.open("{}/NOBS_{}.tif".format(outdir,outfile), 'w', **profile) as dst:
dst.write(S.astype(rasterio.uint16), 1)
#!/usr/bin/env python
# This script computes the Permanent snow area AKA "PSL" in CoSIMS from annnual snow map outputs
# The PSL of a given year is defined as pixels where snow was always observed when possible from May 01 to Sep 01 of this year
# In HAL the dependencies are loaded with module load lis/develop
# Author: Simon Gascoin
import rasterio
import numpy as np
import os,sys
import numpy as np
# input file is input_dates.txt
# Example: f="/work/OT/siaa/Theia/Neige/SNOW_ANNUAL_MAP_LIS_1.5/S2_with_L8_Densification//T32TNS_20170901_20180831/tmpdir/multitemp_cloud_mask.vrt"
f=sys.argv[1]
fdir=os.path.split(os.path.split(f)[0])[0]
fdates=fdir+os.sep+"input_dates.txt"
d=np.loadtxt(fdates)
# Creates May 01 for the year
May01=np.round(d[-1]/10000)*10000+501
# find index of closest date
idx = (np.abs(d - May01)).argmin()
srcObs=rasterio.open(f, 'r')
srcSnow=rasterio.open(os.path.split(f)[0]+os.sep+"multitemp_snow_mask.vrt", 'r')
n=srcObs.meta["count"]
# read only from the date which is the closest to May 01 up to Sep 01
snow = srcSnow.read(range(idx,n))
obs = srcObs.read(range(idx,n))
# pixels in p are either snow (1), cloud (1) or no snow (0)
p=snow+obs
# PSL corresponds to pixels always flagged as snow or cloud
PSL=np.all(p,axis=0)
# output file suffix
outfile=os.path.split(fdir)[1]
# export PSL...tif in the parent folder of the input file
with rasterio.Env():
profile = srcObs.profile
profile.update(
dtype=rasterio.uint8,
driver='GTiff',
count=1)
with rasterio.open("{}/PSL_{}.tif".format(fdir,outfile), 'w', **profile) as dst:
dst.write(PSL.astype(rasterio.uint8), 1)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment