add scripts to compute PSL and NOBS from Theia products

#!/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)