From ca63c140bac1ce0b8214d8551181226909cacc22 Mon Sep 17 00:00:00 2001 From: matheushtavaress Date: Tue, 13 Feb 2024 16:20:12 +0100 Subject: [PATCH] chla by OWTs --- .gitignore | 260 +-- LICENSE | 32 +- README.md | 90 +- getpak/__init__.py | 2 +- getpak/automation.py | 24 +- getpak/commons.py | 29 + getpak/data/means_OWT_Spyrakos_S2A_B1-7.json | 119 ++ getpak/data/means_OWT_Spyrakos_S2A_B2-7.json | 106 ++ getpak/data/s2_proj_ref.json | 264 +-- getpak/inversion_functions.py | 167 +- getpak/raster.py | 1726 +++++++++++------- getpak/validation.py | 210 +-- main.py | 132 +- requirements.txt | 4 +- setup.py | 76 +- 15 files changed, 1962 insertions(+), 1279 deletions(-) create mode 100644 getpak/data/means_OWT_Spyrakos_S2A_B1-7.json create mode 100644 getpak/data/means_OWT_Spyrakos_S2A_B2-7.json diff --git a/.gitignore b/.gitignore index e4a1406..42a8001 100644 --- a/.gitignore +++ b/.gitignore @@ -1,130 +1,130 @@ -.idea -# Byte-compiled / optimized / DLL files -__pycache__/ -*.py[cod] -*$py.class - -# C extensions -*.so - -# Distribution / packaging -.Python -build/ -develop-eggs/ -dist/ -downloads/ -eggs/ -.eggs/ -lib/ -lib64/ -parts/ -sdist/ -var/ -wheels/ -pip-wheel-metadata/ -share/python-wheels/ -*.egg-info/ -.installed.cfg -*.egg -MANIFEST - -# PyInstaller -# Usually these files are written by a python script from a template -# before PyInstaller builds the exe, so as to inject date/other infos into it. -*.manifest -*.spec - -# Installer logs -pip-log.txt -pip-delete-this-directory.txt - -# Unit test / coverage reports -htmlcov/ -.tox/ -.nox/ -.coverage -.coverage.* -.cache -nosetests.xml -coverage.xml -*.cover -*.py,cover -.hypothesis/ -.pytest_cache/ - -# Translations -*.mo -*.pot - -# Django stuff: -*.log -local_settings.py -db.sqlite3 -db.sqlite3-journal - -# Flask stuff: -instance/ -.webassets-cache - -# Scrapy stuff: -.scrapy - -# Sphinx documentation -docs/_build/ - -# PyBuilder -target/ - -# Jupyter Notebook -.ipynb_checkpoints - -# IPython -profile_default/ -ipython_config.py - -# pyenv -.python-version - -# pipenv -# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. -# However, in case of collaboration, if having platform-specific dependencies or dependencies -# having no cross-platform support, pipenv may install dependencies that don't work, or not -# install all needed dependencies. -#Pipfile.lock - -# PEP 582; used by e.g. github.com/David-OConnor/pyflow -__pypackages__/ - -# Celery stuff -celerybeat-schedule -celerybeat.pid - -# SageMath parsed files -*.sage.py - -# Environments -.env -.venv -env/ -venv/ -ENV/ -env.bak/ -venv.bak/ - -# Spyder project settings -.spyderproject -.spyproject - -# Rope project settings -.ropeproject - -# mkdocs documentation -/site - -# mypy -.mypy_cache/ -.dmypy.json -dmypy.json - -# Pyre type checker -.pyre/ +.idea +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +pip-wheel-metadata/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +.python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ diff --git a/LICENSE b/LICENSE index ff06d1a..f368be3 100644 --- a/LICENSE +++ b/LICENSE @@ -1,16 +1,16 @@ -GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 - -Copyright (C) 2023, David Guimaraes. All rights reserved. - -This program is free software: you can redistribute it and/or modify -it under the terms of the GNU 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 General Public License for more details. - -You should have received a copy of the GNU General Public License -along with this program. If not, see . +GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 + +Copyright (C) 2023, David Guimaraes. All rights reserved. + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU 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 General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . diff --git a/README.md b/README.md index e4d67b4..de7d8ec 100644 --- a/README.md +++ b/README.md @@ -1,46 +1,46 @@ -## GET-pak -**G**éosciences **E**nvironnement **T**oulouse - **P**rocessing and **a**nalysis Wor**k**bench - -GETpak aims to provide tools for Sentinel-2, Sentinel-3, GeoTIFF, NetCDF and vector data manipulation and validation. - - - _..._ - .' '. _ - / .-""-\ _/ \ - .-| /:. | | | - | \ |:. /.-'-./ - | .-'-;:__.' =/ ,ad8888ba, 88888888888 888888888888 88 - .'= *=|CNES _.=' d8"' `"8b 88 88 88 - / _. | ; d8' 88 88 88 - ;-.-'| \ | 88 88aaaaa 88 8b,dPPYba, ,adPPYYba, 88 ,d8 - / | \ _\ _\ 88 88888 88""""" 88 aaaaaa 88P' "8a "" `Y8 88 ,a8" - \__/'._;. ==' ==\ Y8, 88 88 88 """""" 88 d8 ,adPPPPP88 8888[ - \ \ | Y8a. .a88 88 88 88b, ,a8" 88, ,88 88`"Yba, - / / / `"Y88888P" 88888888888 88 88`YbbdP"' `"8bbdP"Y8 88 `Y8a - /-._/-._/ 88 - \ `\ \ 88 - `-._/._/ -### Setup -⚠️ GDAL is a requirement for the installation, therefore, -usage of a conda environment -([Anaconda.org](https://www.anaconda.com/products/individual)) -is strongly recommended. Unless you know what you are doing (-: - -## Installation -Create a Conda environment (python versions above 3.9 were not tested but they should also be compatible): -``` -conda create --name getpak python=3.9 -``` -Activate your conda env: -``` -conda activate getpak -``` -With your conda env activated, install GDAL before installing `getpak` to avoid dependecy errors: -``` -conda install -c conda-forge gdal -``` -Clone `getpak` repo to the desired location, enter it and type: -``` -pip install -e . -``` +## GET-pak +**G**éosciences **E**nvironnement **T**oulouse - **P**rocessing and **a**nalysis Wor**k**bench + +GETpak aims to provide tools for Sentinel-2, Sentinel-3, GeoTIFF, NetCDF and vector data manipulation and validation. + + + _..._ + .' '. _ + / .-""-\ _/ \ + .-| /:. | | | + | \ |:. /.-'-./ + | .-'-;:__.' =/ ,ad8888ba, 88888888888 888888888888 88 + .'= *=|CNES _.=' d8"' `"8b 88 88 88 + / _. | ; d8' 88 88 88 + ;-.-'| \ | 88 88aaaaa 88 8b,dPPYba, ,adPPYYba, 88 ,d8 + / | \ _\ _\ 88 88888 88""""" 88 aaaaaa 88P' "8a "" `Y8 88 ,a8" + \__/'._;. ==' ==\ Y8, 88 88 88 """""" 88 d8 ,adPPPPP88 8888[ + \ \ | Y8a. .a88 88 88 88b, ,a8" 88, ,88 88`"Yba, + / / / `"Y88888P" 88888888888 88 88`YbbdP"' `"8bbdP"Y8 88 `Y8a + /-._/-._/ 88 + \ `\ \ 88 + `-._/._/ +### Setup +⚠️ GDAL is a requirement for the installation, therefore, +usage of a conda environment +([Anaconda.org](https://www.anaconda.com/products/individual)) +is strongly recommended. Unless you know what you are doing (-: + +## Installation +Create a Conda environment (python versions above 3.9 were not tested but they should also be compatible): +``` +conda create --name getpak python=3.9 +``` +Activate your conda env: +``` +conda activate getpak +``` +With your conda env activated, install GDAL before installing `getpak` to avoid dependecy errors: +``` +conda install -c conda-forge gdal +``` +Clone `getpak` repo to the desired location, enter it and type: +``` +pip install -e . +``` (Setup in edit mode is strongly recommended until a valid stable release is added to the Python Package Index - PyPI). \ No newline at end of file diff --git a/getpak/__init__.py b/getpak/__init__.py index 8dbfdad..cba8e59 100644 --- a/getpak/__init__.py +++ b/getpak/__init__.py @@ -1 +1 @@ -__version__ = '0.0.3' \ No newline at end of file +__version__ = '0.0.4' \ No newline at end of file diff --git a/getpak/automation.py b/getpak/automation.py index 6eb2c5c..5a2bfac 100644 --- a/getpak/automation.py +++ b/getpak/automation.py @@ -1,12 +1,12 @@ -from getpak.commons import Utils as U -from getpak import raster as R - - -class Pipelines: - def __int__(self): - print(f'mock: starting logger.') # TODO: call logger from utils - - def get_grs_metadict(self, path2grs=None): - mocklist = ['path/to/img1', 'path/to/img2', 'path/to/img3'] - print(path2grs, mocklist) - pass +from getpak.commons import Utils as U +from getpak import raster as R + + +class Pipelines: + def __int__(self): + print(f'mock: starting logger.') # TODO: call logger from utils + + def get_grs_metadict(self, path2grs=None): + mocklist = ['path/to/img1', 'path/to/img2', 'path/to/img3'] + print(path2grs, mocklist) + pass diff --git a/getpak/commons.py b/getpak/commons.py index 10bb9fc..36279e8 100644 --- a/getpak/commons.py +++ b/getpak/commons.py @@ -354,6 +354,7 @@ def sch_date_matchups(fst_dates, snd_dates, fst_tile_list, snd_tile_list): """ Function to search for the matchup dates of two sets of images given two sets of dates. This function also writes the directories of the matchups for each date + Parameters ---------- fst_dates: list of dates of the first set of images @@ -380,6 +381,34 @@ def sch_date_matchups(fst_dates, snd_dates, fst_tile_list, snd_tile_list): return matches, str_matches, dates + @staticmethod + def find_outliers_IQR(values): + """ + Function to remove outliers, using the same methods as in boxplot (interquartile distance) + + Parameters + ---------- + values: numpy array with data + + Returns + ------- + outliers: a numpy array of the removed outliers + clean_array: a numpy array, with the same size, as values, without the outliers + """ + aux = values[np.where(np.isnan(values) == False)] + clean_array = values + q1 = np.quantile(aux, 0.25) + q3 = np.quantile(aux, 0.75) + iqr = q3 - q1 + # finding upper and lower whiskers + upper_bound = q3 + (1.5 * iqr) + lower_bound = q1 - (1.5 * iqr) + # array of outliers and decontaminated array + outliers = aux[(aux <= lower_bound) | (aux >= upper_bound)] + clean_array[(clean_array <= lower_bound) | (clean_array >= upper_bound)] = np.nan + + return outliers, clean_array + class DefaultDicts: clustering_methods = {'M0': ['Oa17_reflectance:float', 'Oa21_reflectance:float'], diff --git a/getpak/data/means_OWT_Spyrakos_S2A_B1-7.json b/getpak/data/means_OWT_Spyrakos_S2A_B1-7.json new file mode 100644 index 0000000..c3aed05 --- /dev/null +++ b/getpak/data/means_OWT_Spyrakos_S2A_B1-7.json @@ -0,0 +1,119 @@ +{ + "OWT1": { + "Band1": 0.029453968950838943, + "Band2": 0.04891433514372519, + "Band3": 0.11212160249608547, + "Band4": 0.05262007066068948, + "Band5": 0.16475407348157756, + "Band6": 0.28721762888201774, + "Band7": 0.30491832038506567 + }, + "OWT2": { + "Band1": 0.10603539172546034, + "Band2": 0.179424806927847, + "Band3": 0.33919755054621287, + "Band4": 0.16082161092333205, + "Band5": 0.13862574513630554, + "Band6": 0.037515174951869124, + "Band7": 0.03837971978897307 + }, + "OWT3": { + "Band1": 0.21874036735321653, + "Band2": 0.29785518211874856, + "Band3": 0.33610354032242856, + "Band4": 0.07193348513062, + "Band5": 0.04445584238305662, + "Band6": 0.015129921933659326, + "Band7": 0.015781660758270352 + }, + "OWT4": { + "Band1": 0.08061231580182746, + "Band2": 0.1464039956977228, + "Band3": 0.29123148749247196, + "Band4": 0.1962171760992044, + "Band5": 0.19033774237342763, + "Band6": 0.047706584395225925, + "Band7": 0.0474906981401198 + }, + "OWT5": { + "Band1": 0.09841662925087848, + "Band2": 0.12323254127221449, + "Band3": 0.17426486448375136, + "Band4": 0.18048626298239914, + "Band5": 0.18519639136218052, + "Band6": 0.12088683444671608, + "Band7": 0.11751647620185995 + }, + "OWT6": { + "Band1": 0.07054506216838155, + "Band2": 0.12839838575883936, + "Band3": 0.2882907200570871, + "Band4": 0.15495633060544461, + "Band5": 0.2240854383561063, + "Band6": 0.06806819055472804, + "Band7": 0.06565587249941311 + }, + "OWT7": { + "Band1": 0.044725391284481164, + "Band2": 0.06937033242424208, + "Band3": 0.16589426397899654, + "Band4": 0.08533063587163703, + "Band5": 0.2748125805582108, + "Band6": 0.17992712111058637, + "Band7": 0.17993967477184608 + }, + "OWT8": { + "Band1": 0.05754043539753355, + "Band2": 0.09978383292023423, + "Band3": 0.2332213086004707, + "Band4": 0.1320458931663692, + "Band5": 0.2614155340186032, + "Band6": 0.10776034575680644, + "Band7": 0.10823265013998272 + }, + "OWT9": { + "Band1": 0.15250889350349348, + "Band2": 0.22494596523738042, + "Band3": 0.34832374334987515, + "Band4": 0.12164382269857796, + "Band5": 0.09001336953881181, + "Band6": 0.02987873682263174, + "Band7": 0.03268546884922942 + }, + "OWT10": { + "Band1": 0.05679849075054226, + "Band2": 0.04834873397096495, + "Band3": 0.10663802798253338, + "Band4": 0.22147851364063356, + "Band5": 0.2721970126950602, + "Band6": 0.1362090091055845, + "Band7": 0.1583302118546812 + }, + "OWT11": { + "Band1": 0.05589755935639065, + "Band2": 0.10393296965516707, + "Band3": 0.21727618518700778, + "Band4": 0.23010875905691786, + "Band5": 0.23994333955948094, + "Band6": 0.07527421490625884, + "Band7": 0.07756697227877678 + }, + "OWT12": { + "Band1": 0.10943482854349587, + "Band2": 0.15010516628485704, + "Band3": 0.2412465387324616, + "Band4": 0.16738084282163648, + "Band5": 0.17877365318091076, + "Band6": 0.07690735356194664, + "Band7": 0.0761516168746915 + }, + "OWT13": { + "Band1": 0.5113425602164394, + "Band2": 0.3440491733493938, + "Band3": 0.10553306599775934, + "Band4": 0.013624524326132722, + "Band5": 0.009064658438149539, + "Band6": 0.006820591538521004, + "Band7": 0.009565426133604274 + } +} \ No newline at end of file diff --git a/getpak/data/means_OWT_Spyrakos_S2A_B2-7.json b/getpak/data/means_OWT_Spyrakos_S2A_B2-7.json new file mode 100644 index 0000000..872395f --- /dev/null +++ b/getpak/data/means_OWT_Spyrakos_S2A_B2-7.json @@ -0,0 +1,106 @@ +{ + "OWT1": { + "Band2": 0.0503987792220929, + "Band3": 0.11552425017377271, + "Band4": 0.05421697578198031, + "Band5": 0.16975400260354293, + "Band6": 0.29593406154217666, + "Band7": 0.31417193067643445 + }, + "OWT2": { + "Band2": 0.2007068347752141, + "Band3": 0.37943062556011625, + "Band4": 0.17989706688023963, + "Band5": 0.15506849359939437, + "Band6": 0.041964944254647814, + "Band7": 0.0429320349303879 + }, + "OWT3": { + "Band2": 0.3812499323812016, + "Band3": 0.43020722724885097, + "Band4": 0.09207372571768593, + "Band5": 0.05690277665114642, + "Band6": 0.019366061295656038, + "Band7": 0.0202002767054591 + }, + "OWT4": { + "Band2": 0.1592407623182449, + "Band3": 0.31676679217914944, + "Band4": 0.21342158424748933, + "Band5": 0.20702663919130837, + "Band6": 0.051889518660272645, + "Band7": 0.05165470340353532 + }, + "OWT5": { + "Band2": 0.13668457656868838, + "Band3": 0.1932875762104568, + "Band4": 0.2001881011097553, + "Band5": 0.20541238599853684, + "Band6": 0.13408281293639185, + "Band7": 0.13034454717617078 + }, + "OWT6": { + "Band2": 0.13814374482574454, + "Band3": 0.31017180965185676, + "Band4": 0.1667174214674157, + "Band5": 0.24109338627958526, + "Band6": 0.07323452464896084, + "Band7": 0.07063911312643695 + }, + "OWT7": { + "Band2": 0.0726182102940209, + "Band3": 0.17366133514431126, + "Band4": 0.08932576569409112, + "Band5": 0.287679142783591, + "Band6": 0.18835120233387123, + "Band7": 0.18836434375011446 + }, + "OWT8": { + "Band2": 0.10587598308509233, + "Band3": 0.24746028090748323, + "Band4": 0.14010775435449768, + "Band5": 0.2773758618799411, + "Band6": 0.11433948978199424, + "Band7": 0.11484062999099141 + }, + "OWT9": { + "Band2": 0.26542575315898925, + "Band3": 0.4110057800958304, + "Band4": 0.1435340403764808, + "Band5": 0.1062115800966023, + "Band6": 0.035255516658044016, + "Band7": 0.03856732961405318 + }, + "OWT10": { + "Band2": 0.051260238132398585, + "Band3": 0.11305964519436512, + "Band4": 0.23481569046350728, + "Band5": 0.28858839815858434, + "Band6": 0.14441135618407916, + "Band7": 0.16786467186706555 + }, + "OWT11": { + "Band2": 0.11008653847384862, + "Band3": 0.23014047611071448, + "Band4": 0.2437328293527651, + "Band5": 0.2541496867605891, + "Band6": 0.07973098221729333, + "Band7": 0.08215948708478941 + }, + "OWT12": { + "Band2": 0.1685504566042737, + "Band3": 0.27089150402985895, + "Band4": 0.1879490105680733, + "Band5": 0.2007417973560874, + "Band6": 0.08635791745164027, + "Band7": 0.08550931399006638 + }, + "OWT13": { + "Band2": 0.7040702654640484, + "Band3": 0.21596533155108158, + "Band4": 0.027881544855159448, + "Band5": 0.018550128781758685, + "Band6": 0.013957817856087538, + "Band7": 0.019574911491864435 + } +} \ No newline at end of file diff --git a/getpak/data/s2_proj_ref.json b/getpak/data/s2_proj_ref.json index a8d40e1..c20f040 100644 --- a/getpak/data/s2_proj_ref.json +++ b/getpak/data/s2_proj_ref.json @@ -1,132 +1,132 @@ -{ - "31TFL": { - "trans": [ - 600000.0, - 20.0, - 0.0, - 5100000.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 31N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",3],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32631\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KLQ": { - "trans": [ - 300000.0, - 20.0, - 0.0, - 7500040.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KMQ": { - "trans": [ - 399960.0, - 20.0, - 0.0, - 7500040.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KNQ": { - "trans": [ - 499980.0, - 20.0, - 0.0, - 7500040.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KNR": { - "trans": [ - 499980.0, - 20.0, - 0.0, - 7600000.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KPQ": { - "trans": [ - 600000.0, - 20.0, - 0.0, - 7500040.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KPR": { - "trans": [ - 600000.0, - 20.0, - 0.0, - 7600000.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - }, - "22KGA": { - "trans": [ - 699960.0, - 20.0, - 0.0, - 7600000.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 22S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-51],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32722\"]]", - "rows": 5490, - "cols": 5490 - }, - "22KHA": { - "trans": [ - 799980.0, - 20.0, - 0.0, - 7600000.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 22S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-51],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32722\"]]", - "rows": 5490, - "cols": 5490 - }, - "23KKQ": { - "trans": [ - 199980.0, - 20.0, - 0.0, - 7500040.0, - 0.0, - -20.0 - ], - "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", - "rows": 5490, - "cols": 5490 - } -} +{ + "31TFL": { + "trans": [ + 600000.0, + 20.0, + 0.0, + 5100000.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 31N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",3],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32631\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KLQ": { + "trans": [ + 300000.0, + 20.0, + 0.0, + 7500040.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KMQ": { + "trans": [ + 399960.0, + 20.0, + 0.0, + 7500040.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KNQ": { + "trans": [ + 499980.0, + 20.0, + 0.0, + 7500040.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KNR": { + "trans": [ + 499980.0, + 20.0, + 0.0, + 7600000.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KPQ": { + "trans": [ + 600000.0, + 20.0, + 0.0, + 7500040.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KPR": { + "trans": [ + 600000.0, + 20.0, + 0.0, + 7600000.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + }, + "22KGA": { + "trans": [ + 699960.0, + 20.0, + 0.0, + 7600000.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 22S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-51],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32722\"]]", + "rows": 5490, + "cols": 5490 + }, + "22KHA": { + "trans": [ + 799980.0, + 20.0, + 0.0, + 7600000.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 22S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-51],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32722\"]]", + "rows": 5490, + "cols": 5490 + }, + "23KKQ": { + "trans": [ + 199980.0, + 20.0, + 0.0, + 7500040.0, + 0.0, + -20.0 + ], + "proj": "PROJCS[\"WGS 84 / UTM zone 23S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-45],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32723\"]]", + "rows": 5490, + "cols": 5490 + } +} diff --git a/getpak/inversion_functions.py b/getpak/inversion_functions.py index d76c7ef..d1784cf 100644 --- a/getpak/inversion_functions.py +++ b/getpak/inversion_functions.py @@ -1,19 +1,34 @@ # This module contains the functions that will be used to invert from reflectances to water quality parameters +# It was imported from the WaterQuality package (https://github.com/cordmaur/WaterQuality) # The dicionary `functions`` will be used by the waterquality module: # function = { # 'name_of_the_parameter': {'function': any_function}, 'units': 'units to be displayed in the report', # 'name_of_the_parameter2': ..... # } -# Any bands can be used to compute the final value. The name of the band must match the internal name used by WaterDetect -# It is enough to put the band name as an argument in the function. +# Any bands can be used to compute the final value # Available bands in Sentinel2 are: # Blue, Green, Red, RedEdg1, RedEdg2, RedEdg3, Nir, Nir2, Mir, Mir2 import numpy as np +#### CDOM #### +# Brezonik et al. 2005 +def cdom_brezonik(Blue, RedEdg2): + cdom = np.exp(1.872 - 0.830 * np.log(Blue/RedEdg2)) + return cdom + +#### Chlorophyll-a #### +# 2-band ratio by Gilerson et al. (2010) +def chl_gilerson2(Red, RedEdg1, a=0.022, b=1.124): + chl = (0.7864 * (RedEdg1 / Red) / a - 0.4245 / a)**b + return chl + +# 3-band ratio by Gilerson et al. (2010) +def chl_gilerson3(Red, RedEdg1, RedEdg2, a=113.36, b=-16.45, c=1.124): + chl = (a * (RedEdg2 *(1 / Red - 1/RedEdg1)) + b)**c + return chl -#### Chlorophyll-a # Gitelson def chl_gitelson(Red, RedEdg1, RedEdg2): chl = 23.1 + 117.4 * (1 / Red - 1 / RedEdg1) * RedEdg2 @@ -24,21 +39,15 @@ def chl_gitelson2(Red, RedEdg1, a=61.324, b=-37.94): chl = a * (RedEdg1 / Red) + b return chl -# OC2 -def chl_OC2(Blue, Green, a=0.2389, b=-1.9369, c=1.7627, d=-3.0777, e=-0.1054): - X = np.log10(Blue / Green) - chl = 10**(a + b*X + c*X**2 + d*X**3 + e*X**4) - return chl - -# 2-band ratio by Gilerson et al. (2010) -def chl_gilerson2(Red, RedEdg1, a=35.75, b=1.124): - chl = (a * (RedEdg1 / Red) - 19.30)**b - return chl - # 2-band semi-analytical by Gons et al. (2003, 2005) -def chl_gons(Red, RedEdg1, RedEdg3, a=1.063, b=0.016, aw664=0.40, aw708=0.70): +def chl_gons(Red, RedEdg1, RedEdg3, a=1.063, b=0.016, aw665=0.40, aw708=0.70): bb = 1.61 * RedEdg3 / (0.082 - 0.6 * RedEdg3) - chl = ((RedEdg1 / Red) * (aw708 + bb) - aw664 - bb**a)/b + chl = ((RedEdg1 / Red) * (aw708 + bb) - aw665 - bb**a)/b + return chl + +# 2-band squared band ration by Gurlin et al. (2011) +def chl_gurlin(Red, RedEdg1, a=25.28, b=14.85, c=-15.18): + chl = (a * (RedEdg1 / Red)**2 + b * (RedEdg1 / Red) + c) return chl # JM Hybride 1 @@ -57,7 +66,19 @@ def chl_h2(Red, RedEdge1, RedEdge2): chl[res >= 0] = 115.794 * RedEdge2[res >= 0] * (1/Red[res >= 0] - 1/RedEdge1[res >= 0]) + 20.678 return chl -#### SPM +# NDCI, Mishra and Mishra (2012) +def chl_ndci(Red, RedEdg1, a=14.039, b=86.115, c=194.325): + index = (RedEdg1 - Red) / (RedEdg1 + Red) + chl = (a + b * index + c * (index*index)) + return chl + +# OC2 +def chl_OC2(Blue, Green, a=0.2389, b=-1.9369, c=1.7627, d=-3.0777, e=-0.1054): + X = np.log10(Blue / Green) + chl = 10**(a + b*X + c*X**2 + d*X**3 + e*X**4) + return chl + +#### SPM #### # Nechad et al. (2010) def nechad(Red, a=610.94, c=0.2324): spm = a * Red / (1 - (Red / c)) @@ -92,7 +113,7 @@ def spm_s3(Red, Nir2, cutoff_value=0.027, cutoff_delta=0.007, low_params=None, h spm = (1 - transition_coef) * low + transition_coef * high return spm -#### Turbidity (FNU) +#### Turbidity #### # Dogliotti def turb_dogliotti(Red, Nir2): """Switching semi-analytical-algorithm computes turbidity from red and NIR band @@ -104,7 +125,7 @@ def turb_dogliotti(Red, Nir2): :return: turbidity in FNU """ - limit_inf, limit_sup = 0.05, 0.07 + limit_inf, limit_sup = 0.05/np.pi, 0.07/np.pi a_low, c_low = 228.1, 0.1641 a_high, c_high = 3078.9, 0.2112 @@ -115,36 +136,51 @@ def turb_dogliotti(Red, Nir2): t_low[Red >= limit_sup] = t_high[Red >= limit_sup] t_low[(Red >= limit_inf) & (Red < limit_sup)] = t_mixing[(Red >= limit_inf) & (Red < limit_sup)] - t_low[t_low > 4000] = 0 + #t_low[t_low > 4000] = 0 return t_low -#### CDOM -# Brezonik et al. 2005 -def cdom_brezonik(Blue, RedEdg2): +def turb_dogliotti_S2(Red, Nir2): + """Switching semi-analytical-algorithm computes turbidity from red and NIR band + following Dogliotti et al., 2015 + The coefficients were recalibrated for Sentinel-2 by Nechad et al. (2016) + :return: turbidity in FNU + """ - cdom = np.exp(1.872 - 0.830 * np.log(Blue/RedEdg2)) + limit_inf, limit_sup = 0.05/np.pi, 0.07/np.pi + a_low, c_low = 610.94, 0.2324 + a_high, c_high = 3030.32, 0.2115 - return cdom + t_low = nechad(Red, a_low, c_low) + t_high = nechad(Nir2, a_high, c_high) + w = (Red - limit_inf) / (limit_sup - limit_inf) + t_mixing = (1 - w) * t_low + w * t_high + t_low[Red >= limit_sup] = t_high[Red >= limit_sup] + t_low[(Red >= limit_inf) & (Red < limit_sup)] = t_mixing[(Red >= limit_inf) & (Red < limit_sup)] + #t_low[t_low > 4000] = 0 + return t_low -functions = { - 'SPM_Nechad': { - 'function': nechad, - 'units': 'mg/l' - }, - - 'SPM_S3': { - 'function': spm_s3, - 'units': 'mg/l' - }, +def turb_conde(Red, a=2.45, b=22.3): + """ + Following the calibration at reservoirs in the Paranapanema river in Condé et al. (2019) + Values between 1.9 and 48 NTU + """ + turb = a * np.exp(b * Red * np.pi) # from rho_w to Rrs + return turb - 'TURB_Dogliotti': { - 'function': turb_dogliotti, - 'units': 'FNU' - }, +# Secchi disk depth +# def secchi_lee(): +# """ +# +# """ +# +# return secchi - 'CHL_Gitelson': { + +functions = { + + 'CHL_Gitelson2': { 'function': chl_gitelson2, 'units': 'mg/m³' }, @@ -154,14 +190,24 @@ def cdom_brezonik(Blue, RedEdg2): 'units': 'mg/m³' }, - 'CHL_Gilerson': { + 'CHL_Gilerson2': { 'function': chl_gilerson2, 'units': 'mg/m³' }, + 'CHL_Gilerson3': { + 'function': chl_gilerson3, + 'units': 'mg/m³' + }, + 'CHL_Gons': { - 'function': chl_gons, - 'units': 'mg/m³' + 'function': chl_gons, + 'units': 'mg/m³' + }, + + 'CHL_Gurlin': { + 'function': chl_gurlin, + 'units': 'mg/m³' }, 'CHL_Hybrid1': { @@ -174,8 +220,43 @@ def cdom_brezonik(Blue, RedEdg2): 'units': 'mg/m³' }, + 'CHL_NDCI': { + 'function': chl_ndci, + 'units': 'mg/m³' + }, + 'CDOM_Brezonik': { 'function': cdom_brezonik, 'units': '', + }, + + # 'SDD_Lee': { + # 'function': secchi_lee, + # 'units': 'm' + # }, + + 'SPM_Nechad': { + 'function': nechad, + 'units': 'mg/l' + }, + + 'SPM_S3': { + 'function': spm_s3, + 'units': 'mg/l' + }, + + 'TURB_Dogliotti': { + 'function': turb_dogliotti, + 'units': 'FNU' + }, + + 'TURB_Dogliotti_S2': { + 'function': turb_dogliotti_S2, + 'units': 'FNU' + }, + + 'TURB_Conde': { + 'function': turb_conde, + 'units': 'NTU' } } \ No newline at end of file diff --git a/getpak/raster.py b/getpak/raster.py index f0fcf6e..07b12ce 100644 --- a/getpak/raster.py +++ b/getpak/raster.py @@ -1,690 +1,1036 @@ -import os -import json -import fiona -import numpy as np -import rasterio -import rasterio.mask -import scipy.ndimage -import importlib_resources -import xarray as xr -from dask.distributed import Client - -from datetime import datetime -from rasterstats import zonal_stats -from rasterio.warp import calculate_default_transform, reproject, Resampling - -from getpak.commons import Utils as u - -class Raster: - """ - Generic class containing methods for matricial manipulations - - Methods - ------- - array2tiff(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS') - Given an input ndarray and the desired projection parameters, create a raster.tif using rasterio. - - array2tiff_gdal(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS') - Given an input ndarray and the desired projection parameters, create a raster.tif using GDT_Float32. - - reproj(in_raster, out_raster, target_crs='EPSG:4326') - Given an input raster.tif reproject it to reprojected.tif using @target_crs - - shp_stats(tif_file, shp_poly, keep_spatial=True, statistics='count min mean max median std') - Given a single-band GeoTIFF file and a vector.shp return statistics inside the polygon. - - extract_px(rasterio_rast, shapefile, rrs_dict, bands) - Given a dict of Rrs and a polygon, to extract the values of pixels from each band - - sam(self, values, single=False) - Given a a set of pixels, uses the Spectral Angle Mapper to generate angle between the Rrs and the OWTs - - classify_owt_px(self, rrs_dict, bands) - Function to classify the OWT of each pixel - - classify_owt(self, rasterio_rast, shapefiles, rrs_dict, bands, min_px=6) - Function to classify the the OWT of pixels inside a shapefile (or a set of shapefiles) - - """ - def __init__(self, parent_log=None): - if parent_log: - self.log = parent_log - else: - INSTANCE_TIME_TAG = datetime.now().strftime('%Y%m%dT%H%M%S') - logfile = os.path.join(os.getcwd(), 'getpak_raster_' + INSTANCE_TIME_TAG + '.log') - # Import CRS projection information from /data/s2_proj_ref.json - s2projdata = importlib_resources.files(__name__).joinpath('data/s2_proj_ref.json') - with s2projdata.open('rb') as fp: - byte_content = fp.read() - self.s2projgrid = json.loads(byte_content) - - # Import OWT means for S2 MSI from /data/means_OWT_S2A_Spyrakos.json - means_owt = importlib_resources.files(__name__).joinpath('data/means_OWT_S2A_Spyrakos.json') - with means_owt.open('rb') as fp: - byte_content = fp.read() - self.owts = dict(json.loads(byte_content)) - - @staticmethod - def array2tiff(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS'): - """ - Given an input ndarray and the desired projection parameters, create a raster.tif using GDT_Float32. - - Parameters - ---------- - @param ndarray_data: Inform if the index should be saved as array in the output folder - @param str_output_file: string of the path of the file to be written - @param transform: rasterio affine transformation matrix (resolution and "upper left" coordinate) - @param projection: projection CRS - @param no_data: the value for no data - @param compression: type of file compression - - @return: None (If all goes well, array2tiff should pass and generate a file inside @str_output_file) - """ - with rasterio.open(fp=str_output_file, - mode='w', - driver='GTiff', - height=ndarray_data.shape[0], - width=ndarray_data.shape[1], - count=1, - dtype=ndarray_data.dtype, - crs=projection, - transform=transform, - nodata=no_data, - options=['COMPRESS=PACKBITS']) as dst: - dst.write(ndarray_data, 1) - - pass - - @staticmethod - def array2tiff_gdal(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS'): - """ - Given an input ndarray and the desired projection parameters, create a raster.tif using GDT_Float32. - - Parameters - ---------- - @param ndarray_data: Inform if the index should be saved as array in the output folder - @param str_output_file: - @param transform: - @param projection: projection CRS - @param no_data: the value for no data - @param compression: type of file compression - - @return: None (If all goes well, array2tiff should pass and generate a file inside @str_output_file) - """ - # Create file using information from the template - outdriver = gdal.GetDriverByName("GTiff") # http://www.gdal.org/gdal_8h.html - # imgs_out = /work/scratch/guimard/grs2spm/ - [cols, rows] = ndarray_data.shape - # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, - # options=['COMPRESS=PACKBITS'] -> https://gdal.org/drivers/raster/gtiff.html#creation-options - outdata = outdriver.Create(str_output_file, rows, cols, 1, gdal.GDT_Float32, options=[compression]) - # Write the array to the file, which is the original array in this example - outdata.GetRasterBand(1).WriteArray(ndarray_data) - # Set a no data value if required - outdata.GetRasterBand(1).SetNoDataValue(no_data) - # Georeference the image - outdata.SetGeoTransform(transform) - # Write projection information - outdata.SetProjection(projection) - # Close the file https://gdal.org/tutorials/raster_api_tut.html#using-create - outdata = None - pass - - @staticmethod - def reproj(in_raster, out_raster, target_crs='EPSG:4326'): - """ - Given an input raster.tif reproject it to reprojected.tif using @target_crs (default = 'EPSG:4326'). - - Parameters - ---------- - @param in_raster: Inform if the index should be saved as array in the output folder - @param out_raster: - @param target_crs: - - @return: None (If all goes well, reproj should pass and generate a reprojected file inside @out_raster) - """ - # Open the input raster file - with rasterio.open(in_raster) as src: - - # Calculate the transformation parameters - transform, width, height = calculate_default_transform( - src.crs, target_crs, src.width, src.height, *src.bounds) - - # Define the output file metadata - kwargs = src.meta.copy() - kwargs.update({ - 'crs': target_crs, - 'transform': transform, - 'width': width, - 'height': height - }) - - # Create the output raster file - with rasterio.open(out_raster, 'w', **kwargs) as dst: - - # Reproject the data - for i in range(1, src.count + 1): - reproject( - source=rasterio.band(src, i), - destination=rasterio.band(dst, i), - src_transform=src.transform, - src_crs=src.crs, - dst_transform=transform, - dst_crs=target_crs, - resampling=Resampling.nearest) - print(f'Done: {out_raster}') - pass - - @staticmethod - def shp_stats(tif_file, shp_poly, keep_spatial=False, statistics='count min mean max median std'): - """ - Given a single-band GeoTIFF file and a vector.shp return statistics inside the polygon. - - Parameters - ---------- - @param tif_file: path to raster.tif file. - @param shp_poly: path to the polygon.shp file. - @param keep_spatial (bool): - True = include the input shp_poly in the output as GeoJSON - False (default) = get only the mini_raster and statistics - @param statistics: what to extract from the shapes, available values are: - - min, max, mean, count, sum, std, median, majority, - minority, unique, range, nodata, percentile. - https://pythonhosted.org/rasterstats/manual.html#zonal-statistics - - @return: roi_stats (dict) containing the extracted statistics inside the region of interest. - """ - # with fiona.open(shp_poly) as src: - # roi_stats = zonal_stats(src, - # tif_file, - # stats=statistics, - # raster_out=True, - # all_touched=True, - # geojson_out=keep_spatial, - # band=1) - # # Original output comes inside a list containing only the output dict: - # return roi_stats[0] - roi_stats = zonal_stats(shp_poly, - tif_file, - stats=statistics, - raster_out=True, - all_touched=True, - geojson_out=keep_spatial, - band=1) - # Original output comes inside a list containing only the output dict: - return roi_stats[0] - - @staticmethod - def extract_px(rasterio_rast, shapefile, rrs_dict, bands): - """ - Given a dict of Rrs and a polygon, to extract the values of pixels from each band - - Parameters - ---------- - rasterio_rast: a rasterio raster open with rasterio.open - shapefile: a polygon opened as geometry using fiona - rrs_dict: a dict containing the Rrs bands - bands: an array containing the bands of the rrs_dict to be extracted - - Returns - ------- - values: an array of dimension n, where n is the number of bands, containing the values inside the polygon - slice: the slice of the polygon (from the rasterio window) - mask_image: the rasterio mask - """ - # rast = rasterio_rast.read(1) - mask_image, _, window_image = rasterio.mask.raster_geometry_mask(rasterio_rast, [shapefile], crop=True) - slices = window_image.toslices() - values = [] - for band in bands: - # subsetting the xarray dataset - subset_data = rrs_dict[band].isel(x=slices[1], y=slices[0]) - # Extract values where mask_image is False - values.append(subset_data.where(~mask_image).values.flatten()) - - return values, slices, mask_image - - def sam(self, rrs, single=False): - """ - Spectral Angle Mapper for OWT classification for a set of pixels - It calculates the angle between the Rrs of a set of pixels and those of the 13 OWT of inland waters (Spyrakos et al., 2018) - Input values are the values of the pixels from B2 to B7, the dict of the OWTs is already stored - Returns the spectral angle between the Rrs of the pixels and each OWT - To classify pixels individually, set single=True - ---------- - """ - if single: - E = rrs / rrs.sum() - else: - med = np.nanmean(rrs, axis=1) - E = med / med.sum() - # norm of the vector - nE = np.linalg.norm(E) - - # Convert owts values to numpy array for vectorized computations - M = np.array([list(val.values()) for val in self.owts.values()]) - nM = np.linalg.norm(M, axis=1) - - # scalar product - num = np.dot(M, E) - den = nM * nE - - angles = np.arccos(num / den) * 180 / np.pi - - return int(np.argmin(angles) + 1) - - def classify_owt_px(self, rrs_dict, bands=['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7']): - """ - Function to classify the OWT of each pixel - - Parameters - ---------- - rrs_dict: a xarray Dataset containing the Rrs bands - bands: an array containing the bands (2 to 7) of the rrs_dict to be extracted - - Returns - ------- - class_px: an array, with the same size as the input bands, with the pixels classified - """ - # Drop the variables that won't be used in the classification - variables_to_drop = [var for var in rrs_dict.variables if var not in bands] - rrs = rrs_dict.drop_vars(variables_to_drop) - # Find non-NaN values - nzero = np.where(~np.isnan(rrs[bands[0]].values)) - # array of OWT class for each pixel - class_px = np.zeros_like(rrs[bands[0]], dtype='uint8') - # array of angles to limit the loop - angles = np.zeros((len(nzero[0])), dtype='uint8') - # loop over each nonzero value in the rrs_dict - pix = np.zeros((len(nzero[0]), len(bands))) - for i in range(len(bands)): - pix[:, i] = rrs[bands[i]].values[nzero] - for i in range(len(nzero[0])): - angles[i] = self.sam(rrs=pix[i, :], single=True) - class_px[nzero] = angles - - return class_px - - def classify_owt(self, rasterio_rast, shapefiles, rrs_dict, bands=['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7'], min_px=6): - """ - Function to classify the the OWT of pixels inside a shapefile (or a set of shapefiles) - - Parameters - ---------- - rasterio_rast: a rasterio raster open with rasterio.open - shapefiles: a polygon (or set of polygons), usually of waterbodies to be classified, opened as geometry using fiona - rrs_dict: a xarray Dataset containing the Rrs bands - bands: an array containing the bands of the rrs_dict to be extracted - min_px: minimum number of pixels in each polygon to operate the classification - - Returns - ------- - class_spt: an array, with the same size as the input bands, with the classified pixels - class_shp: an array with the same length as the shapefiles, with a OWT class for each polygon - """ - class_spt = np.zeros(rrs_dict[bands[0]].shape, dtype='int32') - class_shp = np.zeros((len(shapefiles)), dtype='int32') - for i, shape in enumerate(shapefiles): - values, slices, mask = self.extract_px(rasterio_rast, shape, rrs_dict, bands) - # Verifying if there are more pixels than the minimum - valid_pixels = np.isnan(values[0]) == False - if np.count_nonzero(valid_pixels) >= min_px: - angle = self.sam(values) - else: - angle = int(0) - - # classifying only the valid pixels inside the polygon - values = np.where(valid_pixels, angle, 0) - # adding to avoid replacing values of cropping by other polygons - class_spt[slices[0], slices[1]] += values.reshape((mask.shape)) - # classification by polygon - class_shp[i] = angle - - return class_spt.astype('uint8'), class_shp.astype('uint8') - - def water_colour(self, rasterio_rast, shapefiles, rrs_dict, bands=['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5'], min_px=6): - """ - Function to calculate the water colour based on the Forel-Ule scale, using the bands of Sentinel-2 MSI, with the - coefficients derived by linear correlation by van der Woerd and Wernand (2018). The different combinations of S2 - bands (10, 20 or 60 m resolution) in the visible spectrum can be used, with the default being at 20 m. - - Parameters - ---------- - rasterio_rast: a rasterio raster open with rasterio.open - shapefiles: a polygon (or set of polygons), usually of waterbodies to be classified, opened as geometry using fiona - rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands - bands: an array containing the bands of the rrs_dict to be extracted, in order from blue to red edge - min_px: minimum number of pixels in each polygon to operate the classification - - Returns - ------- - colour: an array, with the same size as the input bands, with the classified pixels - colour_spt: an array with the same length as the shapefiles, with the mean colour (angle) for each polygon - """ - colour_spt = np.zeros(rrs_dict[bands[0]].shape, dtype='float16') - colour_shp = np.zeros((len(shapefiles)), dtype='float16') - - for i, shape in enumerate(shapefiles): - values, slices, mask = self.extract_px(rasterio_rast, shape, rrs_dict, bands) - # Verifying if there are more pixels than the minimum - valid_pixels = np.isnan(values[0]) == False - if np.count_nonzero(valid_pixels) >= min_px: - med = np.nanmean(values, axis=1) - # calculation of the CIE tristimulus (no intercepts): - if len(bands) == 3: - X = 12.040 * med[0] + 53.696 * med[1] + 32.087 * med[2] - Y = 23.122 * med[0] + 65.702 * med[1] + 16.830 * med[2] - Z = 61.055 * med[0] + 1.778 * med[1] + 0.015 * med[2] - delta = lambda x: -164.83 * (alpha / 100) ** 5 + 1139.90 * (alpha / 100) ** 4 - 3006.04 * ( - alpha / 100) ** 3 + 3677.75 * (alpha / 100) ** 2 - 1979.71 * (alpha / 100) + 371.38 - else if len(bands) == 4: - X = 12.040 * med[0] + 53.696 * med[1] + 32.028 * med[2] + 0.529 * med[3] - Y = 23.122 * med[0] + 65.702 * med[1] + 16.808 * med[2] + 0.192 * med[3] - Z = 61.055 * med[0] + 1.778 * med[1] + 0.015 * med[2] + 0.000 * med[3] - delta = lambda x : -161.23 * (alpha / 100) ** 5 + 1117.08 * (alpha / 100) ** 4 - 2950.14 * ( - alpha / 100) ** 3 + 3612.17 * (alpha / 100) ** 2 - 1943.57 * (alpha / 100) + 364.28 - else if len(bands) == 5: - X = 11.756 * med[0] + 6.423 * med[1] + 53.696 * med[2] + 32.028 * med[3] + 0.529 * med[4] - Y = 1.744 * med[0] + 22.289 * med[1] + 65.702 * med[2] + 16.808 * med[3] + 0.192 * med[4] - Z = 62.696 * med[0] + 31.101 * med[1] + 1.778 * med[2] + 0.015 * med[3] + 0.000 * med[4] - delta = lambda x: -65.74 * (alpha / 100) ** 5 + 477.16 * (alpha / 100) ** 4 - 1279.99 * ( - alpha / 100) ** 3 + 1524.96 * (alpha / 100) ** 2 - 751.59 * (alpha / 100) + 116.56 - else: - print("Error in the number of bands provided") - colour_spt, colour_shp = None, None - # normalisation of the tristimulus in 2 coordinates - x = X / (X + Y + Z) - y = Y / (X + Y + Z) - # hue angle: - alpha = (np.arctan2(y - 1 / 3, x - 1 / 3)) * 180 / np.pi % 360 - # correction for multispectral information - alpha_corrected = alpha + delta(alpha) - cls = int(Raster._Forel_Ule_scale(alpha_corrected)) - else: - cls = int(0) - - # classifying only the valid pixels inside the polygon - values = np.where(valid_pixels, cls, 0) - # addition to avoid replacing values by cropping by other polygons - colour_spt[slices[0], slices[1]] += values.reshape((mask.shape)) - # classification by polygon - colour_shp[i] = cls - - return colour_spt, colour_shp - - @staticmethod - def _Forel_Ule_scale(angle): - """ - Calculates the Forel-Ule water colour scale, depending on the hue angle calculated from the Sentinel-2 MSI data, - based on the classification by Novoa et al. (2013) - Parameters - ---------- - angles: the hue angle, in degrees - - Returns - ------- - The water colour class (1-21) - """ - mapping = [(21.0471, 21), (24.4487, 20), (28.2408, 19), (32.6477, 18), (37.1698, 17), (42.3707, 16), - (47.8847, 15), (53.4431, 14), (59.4234, 13), (64.9378, 12), (70.9617, 11), (78.1648, 10), - (88.5017, 9), (99.5371, 8), (118.5208, 7), (147.4148, 6), (178.7020, 5), (202.8305, 4), - (217.1473, 3), (224.8037, 2)] - score = lambda s: next((L for x, L in mapping if s < x), 1) - - return(score(angle)) - - - def chlorophylla(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0): - """ - Function to calculate the chlorophyll-a concentration (chla) based on the optical water type (OWT) - The functions are the ones recomended by Carrea et al. (2023) and Neil et al. (2019, 2020), and are coded in - inversion_functions.py - - Parameters - ---------- - rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands - class_owt_spt: an array, with the same size as the input bands, with the classified pixels - - Returns - ------- - chla: an array, with the same size as the input bands, with the classified pixels - """ - import getpak.inversion_functions as ifunc - chla = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float16') - - # chla functions for each OWT - classes = [1, 4, 5, 6] - index = np.where(np.isin(class_owt_spt, classes)) - if len(index[0]>0): - chla[index] = ifunc.functions['CHL_Gons']['function'](Red=rrs_dict['Rrs_B4'].values[index], - RedEdg1=rrs_dict['Rrs_B5'].values[index], - RedEdg3=rrs_dict['Rrs_B7'].values[index]) - classes = [2, 7, 8, 11, 12] - index = np.where(np.isin(class_owt_spt, classes)) - if len(index[0] > 0): - chla[index] = ifunc.functions['CHL_Gilerson']['function'](Red=rrs_dict['Rrs_B4'].values[index], - RedEdg1=rrs_dict['Rrs_B5'].values[index]) - classes = [3, 9, 10, 13] - index = np.where(np.isin(class_owt_spt, classes)) - if len(index[0] > 0): - chla[index] = ifunc.functions['CHL_OC2']['function'](Blue=rrs_dict['Rrs_B2'].values[index], - Green=rrs_dict['Rrs_B3'].values[index]) - - # removing espurious values - if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): - chla[chla >= upper_lim] = np.nan - chla[chla <= lower_lim] = np.nan - - return chla - - def turbidity(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0): - """ - Function to calculate the turbidity (turb) based on the optical water type (OWT) - The functions are the ones recomended by Carrea et al. (2023) - Parameters - ---------- - rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands - class_owt_spt: an array, with the same size as the input bands, with the classified pixels - - Returns - ------- - turb: an array, with the same size as the input bands, with the classified pixels - """ - import getpak.inversion_functions as ifunc - turb = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float16') - - # chla functions for each OWT - classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13] - index = np.where(np.isin(class_owt_spt, classes)) - if len(index[0] > 0): - turb[index] = ifunc.functions['TURB_Dogliotti']['function'](Red=rrs_dict['Rrs_B4'].values[index], - Nir2=rrs_dict['Rrs_B8A'].values[index]) - - # removing espurious values - if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): - turb[turb >= upper_lim] = np.nan - turb[turb <= lower_lim] = np.nan - - return turb - - def spm(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0): - """ - Function to calculate the suspended particulate matter (SPM) based on the optical water type (OWT) - The functions are the ones recomended by Carrea et al. (2023) - Parameters - ---------- - rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands - class_owt_spt: an array, with the same size as the input bands, with the classified pixels - - Returns - ------- - spm: an array, with the same size as the input bands, with the classified pixels - """ - import getpak.inversion_functions as ifunc - spm = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float16') - - # chla functions for each OWT - classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13] - index = np.where(np.isin(class_owt_spt, classes)) - if len(index[0] > 0): - spm[index] = ifunc.functions['SPM_S3']['function'](Red=rrs_dict['Rrs_B4'].values[index], - Nir2=rrs_dict['Rrs_B8A'].values[index]) - - # removing espurious values - if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): - spm[spm >= upper_lim] = np.nan - spm[spm <= lower_lim] = np.nan - - return spm - -class GRS: - """ - Core functionalities to handle GRS files - - Methods - ------- - metadata(grs_file_entry) - Given a GRS string element, return file metadata extracted from its name. - """ - - def __init__(self, parent_log=None): - if parent_log: - self.log = parent_log - - @staticmethod - def metadata(grs_file_entry): - """ - Given a GRS file return metadata extracted from its name: - - Parameters - ---------- - @param grs_file_entry (str or pathlike obj): path that leads to the GRS.nc file. - - @return: metadata (dict) containing the extracted info, available keys are: - input_file, basename, mission, prod_lvl, str_date, pydate, year, - month, day, baseline_algo_version, relative_orbit, tile, - product_discriminator, cloud_cover, grs_ver. - - Reference - --------- - Given the following file: - /root/23KMQ/2021/05/21/S2A_MSIL1C_20210521T131241_N0300_R138_T23KMQ_20210521T163353_cc020_v15.nc - - S2A : (MMM) is the mission ID(S2A/S2B) - MSIL1C : (MSIXXX) Product procesing level - 20210521T131241 : (YYYYMMDDTHHMMSS) Sensing start time - N0300 : (Nxxyy) Processing Baseline number - R138 : Relative Orbit number (R001 - R143) - T23KMQ : (Txxxxx) Tile Number - 20210521T163353 : Product Discriminator - cc020 : GRS Cloud cover estimation (0-100%) - v15 : GRS algorithm baseline version - - For GRS version >=2.0, the naming does not include the cloud cover estimation nor the GRS version - - Further reading: - Sentinel-2 MSI naming convention: - URL = https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-2-msi/naming-convention - """ - metadata = {} - basefile = os.path.basename(grs_file_entry) - splt = basefile.split('_') - if len(splt)==9: - mission,proc_level,date_n_time,proc_ver,r_orbit,tile,prod_disc,cc,aux = basefile.split('_') - ver,_ = aux.split('.') - else: - mission,proc_level,date_n_time,proc_ver,r_orbit,tile,aux = basefile.split('_') - prod_disc,_ = aux.split('.') - cc, ver = ['NA', 'v20'] - file_event_date = datetime.strptime(date_n_time, '%Y%m%dT%H%M%S') - yyyy = f"{file_event_date.year:02d}" - mm = f"{file_event_date.month:02d}" - dd = f"{file_event_date.day:02d}" - - metadata['input_file'] = grs_file_entry - metadata['basename'] = basefile - metadata['mission'] = mission - metadata['prod_lvl'] = proc_level - metadata['str_date'] = date_n_time - metadata['pydate'] = file_event_date - metadata['year'] = yyyy - metadata['month'] = mm - metadata['day'] = dd - metadata['baseline_algo_version'] = proc_ver - metadata['relative_orbit'] = r_orbit - metadata['tile'] = tile - metadata['product_discriminator'] = prod_disc - metadata['cloud_cover'] = cc - metadata['grs_ver'] = ver - - return metadata - - @staticmethod - def get_grs_dict(grs_nc_file, grs_version='v15'): - """ - Opens the GRS netCDF files using the xarray library and dask, returning a DataArray containing only the Rrs bands - Parameters - ---------- - grs_nc_file: the path to the GRS file - grs_version: a string with GRS version ('v15' or 'v20' for version 2.0.5) - - Returns - ------- - The xarray DataArray containing the 11 Rrs bands, named as 'Rrs_B*' - """ - # list of bands - bands = ['Rrs_B1', 'Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', - 'Rrs_B7', 'Rrs_B8', 'Rrs_B8A', 'Rrs_B11', 'Rrs_B12'] - if grs_version == 'v15': - ds = xr.open_dataset(grs_nc_file, engine="netcdf4", decode_coords='all', chunks={'y': 610, 'x': 610}) - # List of variables to keep - if 'Rrs_B1' in ds.variables: - variables_to_keep = bands - # Drop the variables you don't want - variables_to_drop = [var for var in ds.variables if var not in variables_to_keep] - grs = ds.drop_vars(variables_to_drop) - elif grs_version == 'v20': - ds = xr.open_dataset(grs_nc_file, chunks={'y': 610, 'x': 610}, engine="netcdf4") - waves = ds['Rrs']['wl'] - subset_dict = {band: ds['Rrs'].sel(wl=waves[i]).drop(['wl', 'x', 'y']) for i, band in enumerate(bands)} - grs = xr.Dataset(subset_dict) - else: - grs = None - ds.close() - - return grs - - def param2tiff(self, ndarray_data, img_ref, output_img, no_data=0, gdal_driver_name="GTiff"): - - # Gather information from the template file - ref_data = gdal.Open(img_ref) - trans = ref_data.GetGeoTransform() - proj = ref_data.GetProjection() - # nodatav = 0 #data.GetNoDataValue() - # Create file using information from the template - outdriver = gdal.GetDriverByName(gdal_driver_name) # http://www.gdal.org/gdal_8h.html - - [cols, rows] = ndarray_data.shape - - print(f'Writing output .tiff') - # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, - # options=['COMPRESS=PACKBITS'] -> https://gdal.org/drivers/raster/gtiff.html#creation-options - outdata = outdriver.Create(output_img, rows, cols, 1, gdal.GDT_Float32, options=['COMPRESS=PACKBITS']) - # Write the array to the file, which is the original array in this example - outdata.GetRasterBand(1).WriteArray(ndarray_data) - # Set a no data value if required - outdata.GetRasterBand(1).SetNoDataValue(no_data) - # Georeference the image - outdata.SetGeoTransform(trans) - # Write projection information - outdata.SetProjection(proj) - - # Closing the files - # https://gdal.org/tutorials/raster_api_tut.html#using-create - data = None - outdata = None - self.log.info('') - pass - -# def proc_grs_wd_intersect(self, grs_dict, ): +import os +import json +from osgeo import gdal +import numpy as np +import rasterio +import rasterio.mask +# import scipy.ndimage +import importlib_resources +import xarray as xr +# from dask.distributed import Client + +from datetime import datetime +from rasterstats import zonal_stats +from rasterio.warp import calculate_default_transform, reproject, Resampling + +from getpak.commons import Utils as u + + +class Raster: + """ + Generic class containing methods for matricial manipulations + + Methods + ------- + array2tiff(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS') + Given an input ndarray and the desired projection parameters, create a raster.tif using rasterio. + + array2tiff_gdal(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS') + Given an input ndarray and the desired projection parameters, create a raster.tif using GDT_Float32. + + reproj(in_raster, out_raster, target_crs='EPSG:4326') + Given an input raster.tif reproject it to reprojected.tif using @target_crs + + shp_stats(tif_file, shp_poly, keep_spatial=True, statistics='count min mean max median std') + Given a single-band GeoTIFF file and a vector.shp return statistics inside the polygon. + + extract_px(rasterio_rast, shapefile, rrs_dict, bands) + Given a dict of Rrs and a polygon, to extract the values of pixels from each band + + sam(self, values, single=False) + Given a a set of pixels, uses the Spectral Angle Mapper to generate angle between the Rrs and the OWTs + + classify_owt_px(self, rrs_dict, bands) + Function to classify the OWT of each pixel + + classify_owt(self, rasterio_rast, shapefiles, rrs_dict, bands, min_px=6) + Function to classify the the OWT of pixels inside a shapefile (or a set of shapefiles) + + """ + + def __init__(self, parent_log=None): + if parent_log: + self.log = parent_log + else: + INSTANCE_TIME_TAG = datetime.now().strftime('%Y%m%dT%H%M%S') + logfile = os.path.join(os.getcwd(), 'getpak_raster_' + INSTANCE_TIME_TAG + '.log') + # Import CRS projection information from /data/s2_proj_ref.json + s2projdata = importlib_resources.files(__name__).joinpath('data/s2_proj_ref.json') + with s2projdata.open('rb') as fp: + byte_content = fp.read() + self.s2projgrid = json.loads(byte_content) + + # Import OWT means for S2 MSI from /data/means_OWT_Spyrakos_S2A_B2-7.json + means_owt = importlib_resources.files(__name__).joinpath('data/means_OWT_Spyrakos_S2A_B1-7.json') + with means_owt.open('rb') as fp: + byte_content = fp.read() + self.owts_B1_7 = dict(json.loads(byte_content)) + + means_owt = importlib_resources.files(__name__).joinpath('data/means_OWT_Spyrakos_S2A_B2-7.json') + with means_owt.open('rb') as fp: + byte_content = fp.read() + self.owts_B2_7 = dict(json.loads(byte_content)) + + # raster + + @staticmethod + def array2tiff(ndarray_data, str_output_file, transform, projection, no_data=-1, compression='COMPRESS=PACKBITS'): + """ + Given an input ndarray and the desired projection parameters, create a raster.tif using GDT_Float32. + + Parameters + ---------- + @param ndarray_data: Inform if the index should be saved as array in the output folder + @param str_output_file: string of the path of the file to be written + @param transform: rasterio affine transformation matrix (resolution and "upper left" coordinate) + @param projection: projection CRS + @param no_data: the value for no data + @param compression: type of file compression + + @return: None (If all goes well, array2tiff should pass and generate a file inside @str_output_file) + """ + with rasterio.open(fp=str_output_file, + mode='w', + driver='GTiff', + height=ndarray_data.shape[0], + width=ndarray_data.shape[1], + count=1, + dtype=ndarray_data.dtype, + crs=projection, + transform=transform, + nodata=no_data, + options=[compression]) as dst: + dst.write(ndarray_data, 1) + + pass + + @staticmethod + def array2tiff_gdal(ndarray_data, str_output_file, transform, projection, no_data=-1, + compression='COMPRESS=PACKBITS'): + """ + Given an input ndarray and the desired projection parameters, create a raster.tif using GDT_Float32. + + Parameters + ---------- + @param ndarray_data: Inform if the index should be saved as array in the output folder + @param str_output_file: + @param transform: + @param projection: projection CRS + @param no_data: the value for no data + @param compression: type of file compression + + @return: None (If all goes well, array2tiff should pass and generate a file inside @str_output_file) + """ + # Create file using information from the template + outdriver = gdal.GetDriverByName("GTiff") # http://www.gdal.org/gdal_8h.html + # imgs_out = /work/scratch/guimard/grs2spm/ + [cols, rows] = ndarray_data.shape + # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, + # options=['COMPRESS=PACKBITS'] -> https://gdal.org/drivers/raster/gtiff.html#creation-options + outdata = outdriver.Create(str_output_file, rows, cols, 1, gdal.GDT_Float32, options=[compression]) + # Write the array to the file, which is the original array in this example + outdata.GetRasterBand(1).WriteArray(ndarray_data) + # Set a no data value if required + outdata.GetRasterBand(1).SetNoDataValue(no_data) + # Georeference the image + outdata.SetGeoTransform(transform) + # Write projection information + outdata.SetProjection(projection) + # Close the file https://gdal.org/tutorials/raster_api_tut.html#using-create + outdata = None + pass + + @staticmethod + def reproj(in_raster, out_raster, target_crs='EPSG:4326'): + """ + Given an input raster.tif reproject it to reprojected.tif using @target_crs (default = 'EPSG:4326'). + + Parameters + ---------- + @param in_raster: Inform if the index should be saved as array in the output folder + @param out_raster: + @param target_crs: + + @return: None (If all goes well, reproj should pass and generate a reprojected file inside @out_raster) + """ + # Open the input raster file + with rasterio.open(in_raster) as src: + # Calculate the transformation parameters + transform, width, height = calculate_default_transform( + src.crs, target_crs, src.width, src.height, *src.bounds) + + # Define the output file metadata + kwargs = src.meta.copy() + kwargs.update({ + 'crs': target_crs, + 'transform': transform, + 'width': width, + 'height': height + }) + + # Create the output raster file + with rasterio.open(out_raster, 'w', **kwargs) as dst: + # Reproject the data + for i in range(1, src.count + 1): + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs=target_crs, + resampling=Resampling.nearest) + print(f'Done: {out_raster}') + pass + + @staticmethod + def shp_stats(tif_file, shp_poly, keep_spatial=False, statistics='count min mean max median std'): + """ + Given a single-band GeoTIFF file and a vector.shp return statistics inside the polygon. + + Parameters + ---------- + @param tif_file: path to raster.tif file. + @param shp_poly: path to the polygon.shp file. + @param keep_spatial (bool): + True = include the input shp_poly in the output as GeoJSON + False (default) = get only the mini_raster and statistics + @param statistics: what to extract from the shapes, available values are: + + min, max, mean, count, sum, std, median, majority, + minority, unique, range, nodata, percentile. + https://pythonhosted.org/rasterstats/manual.html#zonal-statistics + + @return: roi_stats (dict) containing the extracted statistics inside the region of interest. + """ + # with fiona.open(shp_poly) as src: + # roi_stats = zonal_stats(src, + # tif_file, + # stats=statistics, + # raster_out=True, + # all_touched=True, + # geojson_out=keep_spatial, + # band=1) + # # Original output comes inside a list containing only the output dict: + # return roi_stats[0] + roi_stats = zonal_stats(shp_poly, + tif_file, + stats=statistics, + raster_out=True, + all_touched=True, + geojson_out=keep_spatial, + band=1) + # Original output comes inside a list containing only the output dict: + return roi_stats[0] + + @staticmethod + def extract_px(rasterio_rast, shapefile, rrs_dict, bands): + """ + Given a dict of Rrs and a polygon, to extract the values of pixels from each band + + Parameters + ---------- + rasterio_rast: a rasterio raster open with rasterio.open + shapefile: a polygon opened as geometry using fiona + rrs_dict: a dict containing the Rrs bands + bands: an array containing the bands of the rrs_dict to be extracted + + Returns + ------- + values: an array of dimension n, where n is the number of bands, containing the values inside the polygon + slice: the slice of the polygon (from the rasterio window) + mask_image: the rasterio mask + """ + # rast = rasterio_rast.read(1) + mask_image, _, window_image = rasterio.mask.raster_geometry_mask(rasterio_rast, [shapefile], crop=True) + slices = window_image.toslices() + values = [] + for band in bands: + # subsetting the xarray dataset + subset_data = rrs_dict[band].isel(x=slices[1], y=slices[0]) + # Extract values where mask_image is False + values.append(subset_data.where(~mask_image).values.flatten()) + + return values, slices, mask_image + + @staticmethod + def extract_function_px(rasterio_rast, shapefiles, data_matrix, fun='median', min_px=6): + """ + Given a numpy matrix of data with the size and projection of a TIFF file opened with rasterio and a polygon, + to extract the values of pixels for each shapefile and returns the values for the desired function + + Parameters + ---------- + rasterio_rast: a rasterio raster open with rasterio.open + shapefiles: set of polygons opened as geometry using fiona + data_matrix: a numpy array with the size and projection of rasterio_rast, with the values to extract + fun: the function to calculate over the data_matrix. Can be one of min, mean, max, median, and std + min_px: minimum number of pixels in each polygon to operate the classification + + Returns + ------- + values_shp: an array with the same length as the shapefiles, with the calculated function for each polygon + """ + + if fun == 'median': + calc = lambda x: np.nanmedian(x) + elif fun == 'mean': + calc = lambda x: np.nanmean(x) + values_shp = np.zeros((len(shapefiles)), dtype='float32') + for i, shape in enumerate(shapefiles): + # extracting the data_matrix by the shapefile + mask_image, _, window_image = rasterio.mask.raster_geometry_mask(rasterio_rast, [shape], crop=True) + slices = window_image.toslices() + values = data_matrix[slices[0], slices[1]][~mask_image] + # Verifying if there are enough pixels to calculate + valid_pixels = np.isnan(values) == False + if np.count_nonzero(valid_pixels) >= min_px: + values_shp[i] = calc(values) + else: + values_shp[i] = np.nan + + return values_shp + + def sam(self, rrs, single=False, mode='B1'): + """ + Spectral Angle Mapper for OWT classification for a set of pixels + It calculates the angle between the Rrs of a set of pixels and those of the 13 OWT of inland waters + (Spyrakos et al., 2018) + Input values are the values of the pixels from B1 or B2 (depends on mode) to B7, the dict of the OWTs is already stored + Returns the spectral angle between the Rrs of the pixels and each OWT + To classify pixels individually, set single=True + ---------- + """ + if single: + E = rrs / rrs.sum() + else: + med = np.nanmedian(rrs, axis=1) + E = med / med.sum() + # norm of the vector + nE = np.linalg.norm(E) + + # Convert OWT values to numpy array for vectorized computations + if mode == 'B1': + M = np.array([list(val.values()) for val in self.owts_B1_7.values()]) + else: + M = np.array([list(val.values()) for val in self.owts_B2_7.values()]) + nM = np.linalg.norm(M, axis=1) + + # scalar product + num = np.dot(M, E) + den = nM * nE + + angles = np.arccos(num / den) + + return angles + + def classify_owt_px(self, rrs_dict, B1=True): + """ + Function to classify the OWT of each pixel + + Parameters + ---------- + rrs_dict: a xarray Dataset containing the Rrs bands + B1: boolean to whether or not use Band 1 when using Sentinel-2 data + + Returns + ------- + class_px: an array, with the same size as the input bands, with the pixels classified with the smallest SAM + angles: a 2-dimensional array, where the first dimension is the number of valid pixels in the rrs_dict, and + the second dimension is the number of OWTs. The values are the spectral angles between the Rrs and the mean Rrs + of each OWT, in each pixel + + """ + if B1: + bands = ['Rrs_B1', 'Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7'] + mode = 'B1' + else: + bands = ['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7'] + mode = 'B2' + # Drop the variables that won't be used in the classification + variables_to_drop = [var for var in rrs_dict.variables if var not in bands + ['x', 'y']] + rrs = rrs_dict.drop_vars(variables_to_drop) + # Find non-NaN values + nzero = np.where(~np.isnan(rrs[bands[-1]].values)) + # array of OWT class for each pixel + class_px = np.zeros_like(rrs[bands[-1]], dtype='uint8') + # array of angles to limit the loop + angles = np.zeros((len(nzero[0]), 13), dtype='float16') + + # creating a new Band 1 by undoing the upsampling of GRS, keeping only the pixels entirely inside water + if B1: + aux = rrs['Rrs_B1'].coarsen(x=3, y=3).mean(skipna=False).interp(x=rrs.x, y=rrs.y, method='nearest').values + rrs['Rrs_B160m'] = (('x', 'y'), aux) + bands = ['Rrs_B160m', 'Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7'] + + # loop over each nonzero value in the rrs_dict + pix = np.zeros((len(nzero[0]), len(bands))) + for i in range(len(bands)): + pix[:, i] = rrs[bands[i]].values[nzero] + for i in range(len(nzero[0])): + angles[i, :] = self.sam(rrs=pix[i, :], mode=mode, single=True) + + # if B1 is being used, there won't be any classification for the pixels without values in B1, so they have to be + # classified using only bands 2 to 7 + if B1: + nodata = np.where(np.isnan(angles[:, 0]))[0] + for i in range(len(nodata)): + angles[nodata[i], :] = self.sam(rrs=pix[nodata[i], 1:], mode='B2', single=True) + class_px[nzero] = np.argmin(angles, axis=1) + 1 + + return class_px, angles + + def classify_owt(self, rasterio_rast, shapefiles, rrs_dict, B1=True, min_px=6): + """ + Function to classify the the OWT of pixels inside a shapefile (or a set of shapefiles) + + Parameters + ---------- + rasterio_rast: a rasterio raster open with rasterio.open + shapefiles: a polygon (or set of polygons), usually of waterbodies to be classified, opened as geometry + using fiona + rrs_dict: a xarray Dataset containing the Rrs bands + B1: boolean to use Band 1 when using Sentinel-2 data + min_px: minimum number of pixels in each polygon to operate the classification + + Returns + ------- + class_spt: an array, with the same size as the input bands, with the classified pixels + class_shp: an array with the same length as the shapefiles, with a OWT class for each polygon + """ + # checking if B1 will be used in the classification + if B1: + bands = ['Rrs_B1', 'Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7'] + mode = 'B1' + else: + bands = ['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', 'Rrs_B7'] + mode = 'B2' + class_spt = np.zeros(rrs_dict[bands[0]].shape, dtype='int32') + class_shp = np.zeros((len(shapefiles)), dtype='int32') + for i, shape in enumerate(shapefiles): + values, slices, mask = self.extract_px(rasterio_rast, shape, rrs_dict, bands) + # Verifying if there are more pixels than the minimum + valid_pixels = np.isnan(values[0]) == False + if np.count_nonzero(valid_pixels) >= min_px: + angle = int(np.argmin(self.sam(values, mode=mode)) + 1) + else: + angle = int(0) + + # classifying only the valid pixels inside the polygon + values = np.where(valid_pixels, angle, 0) + # adding to avoid replacing values of cropping by other polygons + class_spt[slices[0], slices[1]] += values.reshape(mask.shape) + # classification by polygon + class_shp[i] = angle + + return class_spt.astype('uint8'), class_shp.astype('uint8') + + def classify_owt_weights(self, class_px, angles, n=3, remove_cls_one=True): + """ + Function to attribute weights to the n-th most important OWTs, based on the squared cossine of SAM + The weights are used for calculating weighted means of the water quality parameters, in order to smooth the + spatial differences between the pixels, and also to remove possible outliers generated by some models + For more information on this approach, please refer to Moore et al. (2001) and Carrea et al. (2023) + + This function uses the results of classify_owt_px as input data + + Parameters + ---------- + class_px: an array, with the same size as the input bands, with the pixels classified with the smallest SAM + angles: a 2-dimensional array, where the first dimension is the number of valid pixels in the rrs_dict, and + the second dimension is the number of OWTs. The values are the spectral angles between the Rrs and the mean Rrs + of each OWT, in each pixel + n = the number of dominant classes to be used to generate the weights + remove_cls_one: logical, remove the pixels classified as OWT 1? + + Returns + ------- + class_px: an array, with the same size as the input bands, with the pixels classified with the smallest SAM + angles + """ + # creating the variables of weights and classes for each pixel, depending on the desired number of classes to + # be used + owt_weights = np.zeros((n, *class_px.shape), dtype='float32') + owt_classes = np.zeros((n, *class_px.shape), dtype='float32') + + # finding where there is no valid pixels in the reflectance data + nzero = np.where(class_px != 0) + + # Create an array of indices + cos_sqrt = np.cos(angles) ** 2 + indices = np.argsort(-cos_sqrt, axis=1) + weights = np.take_along_axis(cos_sqrt, indices[:, 0:n:], axis=1) + + # Assigning values to the matrices + for i in range(n): + owt_weights[i, nzero[0], nzero[1]] = weights[:, i] + owt_classes[i, nzero[0], nzero[1]] = indices[:, i] + 1 # summing one due to positioning starting in 1 + + # removing OWT 1: + if remove_cls_one == True: + ones = np.where(indices[:, 0] == 0)[0] + owt_weights[:, nzero[0][ones], nzero[1][ones]] = np.nan + + # # removing the zeros to avoid division by 0 + # owt_weights[np.where(owt_weights == 0)] = np.nan + # owt_classes[np.where(owt_classes == 0)] = np.nan + + return owt_classes, owt_weights + + def cdom(self, rrs_dict, class_owt_spt, upper_lim=50, lower_lim=0): + """ + Function to calculate the coloured dissolved organic matter (CDOM) based on the optical water type (OWT) + + Parameters + ---------- + rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands + class_owt_spt: an array, with the same size as the input bands, with the OWT pixels + + Returns + ------- + cdom: an array, with the same size as the input bands, with the modeled values + """ + import getpak.inversion_functions as ifunc + cdom = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float32') + + # create a matrix for the outliers + if not hasattr(self, 'out'): + self.out = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='uint8') + + # spm functions for each OWT + classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + cdom[index] = ifunc.functions['CDOM_Brezonik']['function'](Blue=rrs_dict['Rrs_B2'].values[index], + RedEdg2=rrs_dict['Rrs_B6'].values[index]) + + # removing espurious values + if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): + out = np.where((cdom < lower_lim) | (cdom > upper_lim)) + cdom[out] = np.nan + self.out[out] = 1 + + out = np.where((cdom == 0) | np.isinf(cdom)) + cdom[out] = np.nan + + return cdom + + def chlorophylla(self, rrs_dict, class_owt_spt, limits=True, alg='owt'): + """ + Function to calculate the chlorophyll-a concentration (chla) based on the optical water type (OWT) + The functions are the ones recomended by Carrea et al. (2023) and Neil et al. (2019, 2020), and are coded in + inversion_functions.py + + Parameters + ---------- + rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands + class_owt_spt: an array, with the same size as the input bands, with the OWT pixels + alg: one the the following algorithms available: owt to use the methodology based on OWT, or gons + + Returns + ------- + chla: an array, with the same size as the input bands, with the modeled values + """ + import getpak.inversion_functions as ifunc + chla = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float32') + + # create a matrix for the outliers if it doesn't exist + if not hasattr(self, 'out'): + self.out = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='uint8') + + if alg == 'owt': + # chla functions for each OWT + # for class 1, using the coefficients calibrated by Neil et al. (2020) for the OWT 1 + # classes = [1] + # index = np.where(np.isin(class_owt_spt, classes)) + # if len(index[0] > 0): + # chla[index] = ifunc.functions['CHL_Gurlin']['function'](Red=rrs_dict['Rrs_B4'].values[index], + # RedEdg1=rrs_dict['Rrs_B5'].values[index], + # a=86.09, b=-517.5, c=886.7) + # if limits: + # lims = [10, 1000] + # out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + # chla[index[0][out], index[1][out]] = np.nan + # self.out[index[0][out], index[1][out]] += 1 + + classes = [1, 6, 10] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_Gons']['function'](Red=rrs_dict['Rrs_B4'].values[index], + RedEdg1=rrs_dict['Rrs_B5'].values[index], + RedEdg3=rrs_dict['Rrs_B7'].values[index], + aw665=0.425, aw708=0.704) + if limits: + lims = [1, 250] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + + classes = [2, 4, 5, 12] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_NDCI']['function'](Red=rrs_dict['Rrs_B4'].values[index], + RedEdg1=rrs_dict['Rrs_B5'].values[index]) + if limits: + lims = [5, 250] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + + # for class 2 and 12, when values of NDCI are >20, use Gilerson instead + classes = [2, 12] + conditions = (np.isin(class_owt_spt, classes)) & (chla > 20) + index = np.where(conditions) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_Gilerson2']['function'](Red=rrs_dict['Rrs_B4'].values[index], + RedEdg1=rrs_dict['Rrs_B5'].values[index]) + if limits: + lims = [5, 250] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + + classes = [7, 8, 11] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_Gilerson2']['function'](Red=rrs_dict['Rrs_B4'].values[index], + RedEdg1=rrs_dict['Rrs_B5'].values[index]) + if limits: + lims = [5, 250] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + + # classes = [] + # index = np.where(np.isin(class_owt_spt, classes)) + # if len(index[0] > 0): + # chla[index] = ifunc.functions['CHL_Gilerson3']['function'](Red=rrs_dict['Rrs_B4'].values[index], + # RedEdg1=rrs_dict['Rrs_B5'].values[index], + # RedEdg2=rrs_dict['Rrs_B6'].values[index]) + # if limits: + # lims = [10, 1000] + # out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + # chla[index[0][out], index[1][out]] = np.nan + # self.out[index[0][out], index[1][out]] += 1 + + classes = [3] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_OC2']['function'](Blue=rrs_dict['Rrs_B2'].values[index], + Green=rrs_dict['Rrs_B3'].values[index], + a=0.1098, b=-0.755, c=-14.12, d=-117, e=-17.76) + if limits: + lims = [0.01, 50] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + + classes = [9] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_OC2']['function'](Blue=rrs_dict['Rrs_B2'].values[index], + Green=rrs_dict['Rrs_B3'].values[index], + a=0.0536, b=7.308, c=116.2, d=412.4, e=463.5) + if limits: + lims = [0.01, 50] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + + classes = [13] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + chla[index] = ifunc.functions['CHL_OC2']['function'](Blue=rrs_dict['Rrs_B2'].values[index], + Green=rrs_dict['Rrs_B3'].values[index], + a=-5020, b=2.9e+04, c=-6.1e+04, d=5.749e+04, e=-2.026e+04) + if limits: + lims = [0.01, 50] + out = np.where((chla[index] < lims[0]) | (chla[index] > lims[1])) + chla[index[0][out], index[1][out]] = np.nan + self.out[index[0][out], index[1][out]] += 1 + else: + chla = ifunc.functions['CHL_Gons']['function'](Red=rrs_dict['Rrs_B4'].values, + RedEdg1=rrs_dict['Rrs_B5'].values, + RedEdg3=rrs_dict['Rrs_B7'].values) + if limits: + lims = [1, 250] + out = np.where((chla < lims[0]) | (chla > lims[1])) + chla[out] = np.nan + self.out[out] += 1 + + # removing espurious values and zeros + out = np.where((chla == 0) | np.isinf(chla)) + chla[out] = np.nan + + return chla + + def turbidity(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0, alg='Dogliotti'): + """ + Function to calculate the turbidity (turb) based on the optical water type (OWT) + Parameters + ---------- + rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands + class_owt_spt: an array, with the same size as the input bands, with the OWT pixels + alg: one the the following algorithms available: Dogliotti or Conde + + Returns + ------- + turb: an array, with the same size as the input bands, with the modeled values + """ + import getpak.inversion_functions as ifunc + turb = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float32') + + # create a matrix for the outliers + if not hasattr(self, 'out'): + self.out = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='uint8') + + # turbidity functions for each OWT + classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + if alg == 'Dogliotti': + turb[index] = ifunc.functions['TURB_Dogliotti_S2']['function'](Red=rrs_dict['Rrs_B4'].values[index], + Nir2=rrs_dict['Rrs_B8A'].values[index]) + elif alg == 'Conde': + turb[index] = ifunc.functions['TURB_Conde']['function'](Red=rrs_dict['Rrs_B4'].values[index]) + + # removing espurious values and zeros + if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): + out = np.where((turb < lower_lim) | (turb > upper_lim)) + turb[out] = np.nan + self.out[out] = 1 + + out = np.where((turb == 0) | np.isinf(turb)) + turb[out] = np.nan + + return turb + + def spm(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0, alg='Hybrid'): + """ + Function to calculate the suspended particulate matter (SPM) based on the optical water type (OWT) + + Parameters + ---------- + rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands + class_owt_spt: an array, with the same size as the input bands, with the OWT pixels + alg: one the the following algorithms available: Hybrid or Nechad + + Returns + ------- + spm: an array, with the same size as the input bands, with the modeled values + """ + import getpak.inversion_functions as ifunc + spm = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float32') + + # create a matrix for the outliers + if not hasattr(self, 'out'): + self.out = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='uint8') + + # spm functions for each OWT + classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + if alg == 'Hybrid': + spm[index] = ifunc.functions['SPM_S3']['function'](Red=rrs_dict['Rrs_B4'].values[index], + Nir2=rrs_dict['Rrs_B8A'].values[index]) + elif alg == 'Nechad': + spm[index] = ifunc.functions['SPM_Nechad']['function'](Red=rrs_dict['Rrs_B4'].values[index]) + + # removing espurious values and zeros + if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): + out = np.where((spm < lower_lim) | (spm > upper_lim)) + spm[out] = np.nan + self.out[out] = 1 + + out = np.where((spm == 0) | np.isinf(spm)) + spm[out] = np.nan + + return spm + + def secchi_dd(self, rrs_dict, class_owt_spt, upper_lim=50, lower_lim=0): + """ + Function to calculate the Secchi disk depth (SPM) based on the optical water type (OWT) + + Parameters + ---------- + rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands + class_owt_spt: an array, with the same size as the input bands, with the OWT pixels + + Returns + ------- + secchi: an array, with the same size as the input bands, with the modeled values + """ + import getpak.inversion_functions as ifunc + secchi = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float32') + + # create a matrix for the outliers + if not hasattr(self, 'out'): + self.out = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='uint8') + + # spm functions for each OWT + classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13] + index = np.where(np.isin(class_owt_spt, classes)) + if len(index[0] > 0): + secchi[index] = ifunc.functions['SDD_Lee']['function'](Red=rrs_dict['Rrs_B4'].values[index]) + + # removing espurious values and zeros + if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)): + out = np.where((secchi < lower_lim) | (secchi > upper_lim)) + secchi[out] = np.nan + self.out[out] = 1 + + out = np.where((secchi == 0) | np.isinf(secchi)) + secchi[out] = np.nan + + return secchi + + @staticmethod + def water_colour(rrs_dict, bands=['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5']): + """ + Function to calculate the water colour of each pixel based on the Forel-Ule scale, using the bands of Sentinel-2 + MSI, with coefficients derived by linear correlation by van der Woerd and Wernand (2018). The different + combinations of S2 bands (10, 20 or 60 m resolution) in the visible spectrum can be used, with the default + being at 20 m. + + Parameters + ---------- + rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands + bands: an array containing the bands of the rrs_dict to be extracted, in order from blue to red edge + + Returns + ------- + colour: an array, with the same size as the input bands, with the classified pixels + """ + # Drop the variables that won't be used in the classification + variables_to_drop = [var for var in rrs_dict.variables if var not in bands] + rrs = rrs_dict.drop_vars(variables_to_drop) + # Find non-NaN values + nzero = np.where(~np.isnan(rrs[bands[0]].values)) + # array of classes to limit the loop + classes = np.zeros((len(nzero[0])), dtype='uint8') + # array of colour class for each pixel + colour = np.zeros_like(rrs[bands[0]], dtype='uint8') + + # calculation of the CIE tristimulus (no intercepts): + # B2 to B4 + if len(bands) == 3: + X = (12.040 * rrs[bands[0]].values[nzero] + 53.696 * rrs[bands[1]].values[nzero] + + 32.087 * rrs[bands[2]].values[nzero]) + Y = (23.122 * rrs[bands[0]].values[nzero] + 65.702 * rrs[bands[1]].values[nzero] + + 16.830 * rrs[bands[2]].values[nzero]) + Z = (61.055 * rrs[bands[0]].values[nzero] + 1.778 * rrs[bands[1]].values[nzero] + + 0.015 * rrs[bands[2]].values[nzero]) + delta = lambda d: -164.83 * (alpha / 100) ** 5 + 1139.90 * (alpha / 100) ** 4 - 3006.04 * ( + alpha / 100) ** 3 + 3677.75 * (alpha / 100) ** 2 - 1979.71 * (alpha / 100) + 371.38 + # B2 to B5 + elif len(bands) == 4: + X = (12.040 * rrs[bands[0]].values[nzero] + 53.696 * rrs[bands[1]].values[nzero] + + 32.028 * rrs[bands[2]].values[nzero] + 0.529 * rrs[bands[3]].values[nzero]) + Y = (23.122 * rrs[bands[0]].values[nzero] + 65.702 * rrs[bands[1]].values[nzero] + + 16.808 * rrs[bands[2]].values[nzero] + 0.192 * rrs[bands[3]].values[nzero]) + Z = (61.055 * rrs[bands[0]].values[nzero] + 1.778 * rrs[bands[1]].values[nzero] + + 0.015 * rrs[bands[2]].values[nzero] + 0.000 * rrs[bands[3]].values[nzero]) + delta = lambda d: -161.23 * (alpha / 100) ** 5 + 1117.08 * (alpha / 100) ** 4 - 2950.14 * ( + alpha / 100) ** 3 + 3612.17 * (alpha / 100) ** 2 - 1943.57 * (alpha / 100) + 364.28 + # B1 to B5 + elif len(bands) == 5: + X = (11.756 * rrs[bands[0]].values[nzero] + 6.423 * rrs[bands[1]].values[nzero] + + 53.696 * rrs[bands[2]].values[nzero] + 32.028 * rrs[bands[3]].values[nzero] + + 0.529 * rrs[bands[4]].values[nzero]) + Y = (1.744 * rrs[bands[0]].values[nzero] + 22.289 * rrs[bands[1]].values[nzero] + + 65.702 * rrs[bands[2]].values[nzero] + 16.808 * rrs[bands[3]].values[nzero] + + 0.192 * rrs[bands[4]].values[nzero]) + Z = (62.696 * rrs[bands[0]].values[nzero] + 31.101 * rrs[bands[1]].values[nzero] + + 1.778 * rrs[bands[2]].values[nzero] + 0.015 * rrs[bands[3]].values[nzero] + + 0.000 * rrs[bands[4]].values[nzero]) + delta = lambda d: -65.74 * (alpha / 100) ** 5 + 477.16 * (alpha / 100) ** 4 - 1279.99 * ( + alpha / 100) ** 3 + 1524.96 * (alpha / 100) ** 2 - 751.59 * (alpha / 100) + 116.56 + else: + print("Error in the number of bands provided") + colour = None + # normalisation of the tristimulus in 2 coordinates + x = X / (X + Y + Z) + y = Y / (X + Y + Z) + # hue angle: + alpha = (np.arctan2(y - 1 / 3, x - 1 / 3)) * 180 / np.pi % 360 + # correction for multispectral information + alpha_corrected = alpha + delta(alpha) + for i in range(len(nzero[0])): + classes[i] = int(Raster._Forel_Ule_scale(alpha_corrected[i])) + # adding the classes to the matrix + colour[nzero] = classes + + return colour + + @staticmethod + def _Forel_Ule_scale(angle): + """ + Calculates the Forel-Ule water colour scale, depending on the hue angle calculated from the Sentinel-2 MSI data, + based on the classification by Novoa et al. (2013) + Parameters + ---------- + angle: the hue angle, in degrees + + Returns + ------- + The water colour class (1-21) + """ + mapping = [(21.0471, 21), (24.4487, 20), (28.2408, 19), (32.6477, 18), (37.1698, 17), (42.3707, 16), + (47.8847, 15), (53.4431, 14), (59.4234, 13), (64.9378, 12), (70.9617, 11), (78.1648, 10), + (88.5017, 9), (99.5371, 8), (118.5208, 7), (147.4148, 6), (178.7020, 5), (202.8305, 4), + (217.1473, 3), (224.8037, 2)] + score = lambda s: next((L for x, L in mapping if s < x), 1) + + return score(angle) + +class GRS: + """ + Core functionalities to handle GRS files + + Methods + ------- + metadata(grs_file_entry) + Given a GRS string element, return file metadata extracted from its name. + """ + + def __init__(self, parent_log=None): + if parent_log: + self.log = parent_log + + @staticmethod + def metadata(grs_file_entry): + """ + Given a GRS file return metadata extracted from its name: + + Parameters + ---------- + @param grs_file_entry: str or pathlike obj that leads to the GRS.nc file. + + @return: metadata (dict) containing the extracted info, available keys are: + input_file, basename, mission, prod_lvl, str_date, pydate, year, + month, day, baseline_algo_version, relative_orbit, tile, + product_discriminator, cloud_cover, grs_ver. + + Reference + --------- + Given the following file: + /root/23KMQ/2021/05/21/S2A_MSIL1C_20210521T131241_N0300_R138_T23KMQ_20210521T163353_cc020_v15.nc + + S2A : (MMM) is the mission ID(S2A/S2B) + MSIL1C : (MSIXXX) Product procesing level + 20210521T131241 : (YYYYMMDDTHHMMSS) Sensing start time + N0300 : (Nxxyy) Processing Baseline number + R138 : Relative Orbit number (R001 - R143) + T23KMQ : (Txxxxx) Tile Number + 20210521T163353 : Product Discriminator + cc020 : GRS Cloud cover estimation (0-100%) + v15 : GRS algorithm baseline version + + For GRS version >=2.0, the naming does not include the cloud cover estimation nor the GRS version + + Further reading: + Sentinel-2 MSI naming convention: + URL = https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-2-msi/naming-convention + """ + metadata = {} + basefile = os.path.basename(grs_file_entry) + splt = basefile.split('_') + if len(splt) == 9: + mission, proc_level, date_n_time, proc_ver, r_orbit, tile, prod_disc, cc, aux = basefile.split('_') + ver, _ = aux.split('.') + else: + mission, proc_level, date_n_time, proc_ver, r_orbit, tile, aux = basefile.split('_') + prod_disc, _ = aux.split('.') + cc, ver = ['NA', 'v20'] + file_event_date = datetime.strptime(date_n_time, '%Y%m%dT%H%M%S') + yyyy = f"{file_event_date.year:02d}" + mm = f"{file_event_date.month:02d}" + dd = f"{file_event_date.day:02d}" + + metadata['input_file'] = grs_file_entry + metadata['basename'] = basefile + metadata['mission'] = mission + metadata['prod_lvl'] = proc_level + metadata['str_date'] = date_n_time + metadata['pydate'] = file_event_date + metadata['year'] = yyyy + metadata['month'] = mm + metadata['day'] = dd + metadata['baseline_algo_version'] = proc_ver + metadata['relative_orbit'] = r_orbit + metadata['tile'] = tile + metadata['product_discriminator'] = prod_disc + metadata['cloud_cover'] = cc + metadata['grs_ver'] = ver + + return metadata + + @staticmethod + def get_grs_dict(grs_nc_file, grs_version='v15'): + """ + Opens the GRS netCDF files using the xarray library and dask, returning a DataArray containing only the + Rrs bands + Parameters + ---------- + grs_nc_file: the path to the GRS file + grs_version: a string with GRS version ('v15' or 'v20' for version 2.0.5) + + Returns + ------- + The xarray DataArray containing the 11 Rrs bands, named as 'Rrs_B*' + """ + # list of bands + bands = ['Rrs_B1', 'Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6', + 'Rrs_B7', 'Rrs_B8', 'Rrs_B8A', 'Rrs_B11', 'Rrs_B12'] + if grs_version == 'v15': + ds = xr.open_dataset(grs_nc_file, engine="netcdf4", decode_coords='all', chunks={'y': 610, 'x': 610}) + # List of variables to keep + if 'Rrs_B1' in ds.variables: + variables_to_keep = bands + # Drop the variables you don't want + variables_to_drop = [var for var in ds.variables if var not in variables_to_keep] + grs = ds.drop_vars(variables_to_drop) + elif grs_version == 'v20': + ds = xr.open_dataset(grs_nc_file, chunks={'y': -1, 'x': -1}, engine="netcdf4") + waves = ds['Rrs']['wl'] + subset_dict = {band: ds['Rrs'].sel(wl=waves[i]).drop(['wl']) for i, band in enumerate(bands)} + grs = xr.Dataset(subset_dict) + else: + grs = None + ds.close() + + return grs + + def param2tiff(self, ndarray_data, img_ref, output_img, no_data=0, gdal_driver_name="GTiff"): + + # Gather information from the template file + ref_data = gdal.Open(img_ref) + trans = ref_data.GetGeoTransform() + proj = ref_data.GetProjection() + # nodatav = 0 #data.GetNoDataValue() + # Create file using information from the template + outdriver = gdal.GetDriverByName(gdal_driver_name) # http://www.gdal.org/gdal_8h.html + + [cols, rows] = ndarray_data.shape + + print(f'Writing output .tiff') + # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, + # options=['COMPRESS=PACKBITS'] -> https://gdal.org/drivers/raster/gtiff.html#creation-options + outdata = outdriver.Create(output_img, rows, cols, 1, gdal.GDT_Float32, options=['COMPRESS=PACKBITS']) + # Write the array to the file, which is the original array in this example + outdata.GetRasterBand(1).WriteArray(ndarray_data) + # Set a no data value if required + outdata.GetRasterBand(1).SetNoDataValue(no_data) + # Georeference the image + outdata.SetGeoTransform(trans) + # Write projection information + outdata.SetProjection(proj) + + # Closing the files + # https://gdal.org/tutorials/raster_api_tut.html#using-create + # data = None + outdata = None + self.log.info('') + pass + +# def proc_grs_wd_intersect(self, grs_dict, ): diff --git a/getpak/validation.py b/getpak/validation.py index 7a8286c..5004058 100644 --- a/getpak/validation.py +++ b/getpak/validation.py @@ -1,105 +1,105 @@ -import sys -import numpy as np -import pandas as pd -from sklearn.metrics import r2_score -from sklearn.metrics import mean_squared_error -from sklearn.metrics import mean_squared_log_error -from sklearn.metrics import mean_absolute_percentage_error - - -class Validation: - - def __init__(self, parent_log=None): - if parent_log: - self.log = parent_log - - @staticmethod - def r2(y_true, y_pred): #1 - return r2_score(y_true, y_pred) - - @staticmethod - def bias(y_true, y_pred): #2 - return np.mean(np.asarray(y_pred) - np.asarray(y_true)) - - @staticmethod - def rmse(y_true, y_pred): #3 - # return np.sqrt(((y_pred - y_true) ** 2).mean()) - return np.sqrt(mean_squared_error( y_true, y_pred )) - - @staticmethod - def rrmse(y_true, y_pred): #4 - # https://www.analyticsvidhya.com/blog/2021/10/evaluation-metric-for-regression-models/ - """ - Model accuracy is: - RRMSE < 10% (Excellent) - RRMSE is between 10% and 20% (Good) - RRMSE is between 20% and 30% (Fair) - RRMSE > 30% (Poor) - """ - num = np.sum(np.square(np.asarray(y_true) - np.asarray(y_pred))) - den = np.sum(np.square(y_pred)) - squared_error = num/den - rrmse = np.sqrt(squared_error) - return rrmse - - # implementation of NRMSE with standard deviation - def nrmse(self, y_true, y_pred): #5 - """ - NRMSE is a good measure when you want to compare the models of different dependent variables or - when the dependent variables are modified (log-transformed or standardized). It overcomes the - scale-dependency and eases comparison between models of different scales or even between datasets. - """ - local_rmse = self.rmse(y_true, y_pred) - nrmse = local_rmse / np.std(y_pred) - return nrmse - - @staticmethod - def rmsle(y_true, y_pred): #6 - return np.sqrt(mean_squared_log_error( y_true, y_pred )) - - # TODO: Add WMAPE - - @staticmethod - def mape(y_true, y_pred, drop_zero=True): #7 - if 0 in y_true and drop_zero: - print('Zero values present in Y-truth set. MAPE will attempt to be performed without them.') - try: - df = pd.DataFrame(data={'y_true':y_true,'y_pred':y_pred}) - zero_tot = len(df.loc[df['y_true'] == 0]) - print(f'Removing {zero_tot}/{len(df)} values.') - df = df[df['y_true'] != 0].copy() - return mean_absolute_percentage_error(df['y_true'], df['y_pred']) - except: - e = sys.exc_info()[0] - print(e) - return mean_absolute_percentage_error(y_true, y_pred) - - @staticmethod - def corr(y_true, y_pred, m='spearman'): # 8 #9 #10 - # Test for same DataFrame - if isinstance(y_true, pd.Series) and isinstance(y_true, pd.Series): - # Compute pairwise correlation of columns, excluding NA/null values. - # https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.corr.html - return y_true.corr(y_pred, method=m) - else: - # Test for other dtypes and convert to DF - try: - df = pd.DataFrame(data={'y_true': y_true, 'y_pred': y_pred}) - return df['y_true'].corr(df['y_pred'], method=m) - except: - e = sys.exc_info()[0] - print(e) - - def get_stats(self, y_true, y_pred): - res = {} - res['R2'] = round(self.r2(y_true, y_pred), 2) - res['BIAS'] = round(self.bias(y_true, y_pred), 2) - res['MAPE'] = round(self.mape(y_true, y_pred), 2) - res['RMSE'] = round(self.rmse(y_true, y_pred), 2) - res['RRMSE'] = round(self.rrmse(y_true, y_pred), 2) - res['NRMSE'] = round(self.nrmse(y_true, y_pred), 2) - res['RMSLE'] = round(self.rmsle(y_true, y_pred), 2) - res['Spearman'] = round(self.corr(y_true, y_pred), 2) # default m = spearman - res['Pearson'] = round(self.corr(y_true, y_pred, m='pearson'), 2) - res['Kendall'] = round(self.corr(y_true, y_pred, m='kendall'), 2) - return res +import sys +import numpy as np +import pandas as pd +from sklearn.metrics import r2_score +from sklearn.metrics import mean_squared_error +from sklearn.metrics import mean_squared_log_error +from sklearn.metrics import mean_absolute_percentage_error + + +class Validation: + + def __init__(self, parent_log=None): + if parent_log: + self.log = parent_log + + @staticmethod + def r2(y_true, y_pred): #1 + return r2_score(y_true, y_pred) + + @staticmethod + def bias(y_true, y_pred): #2 + return np.mean(np.asarray(y_pred) - np.asarray(y_true)) + + @staticmethod + def rmse(y_true, y_pred): #3 + # return np.sqrt(((y_pred - y_true) ** 2).mean()) + return np.sqrt(mean_squared_error( y_true, y_pred )) + + @staticmethod + def rrmse(y_true, y_pred): #4 + # https://www.analyticsvidhya.com/blog/2021/10/evaluation-metric-for-regression-models/ + """ + Model accuracy is: + RRMSE < 10% (Excellent) + RRMSE is between 10% and 20% (Good) + RRMSE is between 20% and 30% (Fair) + RRMSE > 30% (Poor) + """ + num = np.sum(np.square(np.asarray(y_true) - np.asarray(y_pred))) + den = np.sum(np.square(y_pred)) + squared_error = num/den + rrmse = np.sqrt(squared_error) + return rrmse + + # implementation of NRMSE with standard deviation + def nrmse(self, y_true, y_pred): #5 + """ + NRMSE is a good measure when you want to compare the models of different dependent variables or + when the dependent variables are modified (log-transformed or standardized). It overcomes the + scale-dependency and eases comparison between models of different scales or even between datasets. + """ + local_rmse = self.rmse(y_true, y_pred) + nrmse = local_rmse / np.std(y_pred) + return nrmse + + @staticmethod + def rmsle(y_true, y_pred): #6 + return np.sqrt(mean_squared_log_error( y_true, y_pred )) + + # TODO: Add WMAPE + + @staticmethod + def mape(y_true, y_pred, drop_zero=True): #7 + if 0 in y_true and drop_zero: + print('Zero values present in Y-truth set. MAPE will attempt to be performed without them.') + try: + df = pd.DataFrame(data={'y_true':y_true,'y_pred':y_pred}) + zero_tot = len(df.loc[df['y_true'] == 0]) + print(f'Removing {zero_tot}/{len(df)} values.') + df = df[df['y_true'] != 0].copy() + return mean_absolute_percentage_error(df['y_true'], df['y_pred']) + except: + e = sys.exc_info()[0] + print(e) + return mean_absolute_percentage_error(y_true, y_pred) + + @staticmethod + def corr(y_true, y_pred, m='spearman'): # 8 #9 #10 + # Test for same DataFrame + if isinstance(y_true, pd.Series) and isinstance(y_true, pd.Series): + # Compute pairwise correlation of columns, excluding NA/null values. + # https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.corr.html + return y_true.corr(y_pred, method=m) + else: + # Test for other dtypes and convert to DF + try: + df = pd.DataFrame(data={'y_true': y_true, 'y_pred': y_pred}) + return df['y_true'].corr(df['y_pred'], method=m) + except: + e = sys.exc_info()[0] + print(e) + + def get_stats(self, y_true, y_pred): + res = {} + res['R2'] = round(self.r2(y_true, y_pred), 2) + res['BIAS'] = round(self.bias(y_true, y_pred), 2) + res['MAPE'] = round(self.mape(y_true, y_pred), 2) + res['RMSE'] = round(self.rmse(y_true, y_pred), 2) + res['RRMSE'] = round(self.rrmse(y_true, y_pred), 2) + res['NRMSE'] = round(self.nrmse(y_true, y_pred), 2) + res['RMSLE'] = round(self.rmsle(y_true, y_pred), 2) + res['Spearman'] = round(self.corr(y_true, y_pred), 2) # default m = spearman + res['Pearson'] = round(self.corr(y_true, y_pred, m='pearson'), 2) + res['Kendall'] = round(self.corr(y_true, y_pred, m='kendall'), 2) + return res diff --git a/main.py b/main.py index a30ba5d..542ca03 100644 --- a/main.py +++ b/main.py @@ -1,66 +1,66 @@ -import time -import argparse - -import getpak -from getpak.commons import Utils as U -from getpak import automation as A - - -""" -Author: David Guimaraes - dvdgmf@gmail.com -""" - -# ,-------------, -# | ENTRY POINT | -# '-------------' -if __name__ == '__main__': - # ,--------------, - # | Start timers | - # '--------------' - U.tic() - t1 = time.perf_counter() - - # ,------, - # | LOGO | - # '------' - U.print_logo() - - # ,-----------------, - # | Present options | - # '-----------------' - parser = argparse.ArgumentParser() - parser.add_argument('-v', '--version', help='Displays current package version.', action='store_true') - parser.add_argument('-mg', help="Given a path, return GRS metadata", default='local', type=str) - args = parser.parse_args().__dict__ # Converts the input arguments from Namespace() to dict - - print('User Input Arguments:') - for key in args: - print(f'{key}: {args[key]}') - - # ,----------------------------------------, - # | Automation pipelines class declaration | - # '----------------------------------------' - gpk_pipe = A.Pipelines() - - # ,---------------, - # | Treat options | - # '---------------' - if args['version']: - print(f'GET-Pak version: {getpak.__version__}') - elif args['mg']: - gpk_pipe.get_grs_metadict(args['mg']) - else: - print('No commands supplied, exiting.\n') - - # ,------------------------------, - # | End timers and report to log | - # '------------------------------' - t_hour, t_min, t_sec, _ = U.tac() - t2 = time.perf_counter() - final_msg_1 = f'Finished in {round(t2 - t1, 2)} second(s).' - final_msg_2 = f'Elapsed execution time: {t_hour}h : {t_min}m : {t_sec}s' - print(final_msg_1) - print(final_msg_2) - # ,-----, - # | END | - # '-----' +import time +import argparse + +import getpak +from getpak.commons import Utils as U +from getpak import automation as A + + +""" +Author: David Guimaraes - dvdgmf@gmail.com +""" + +# ,-------------, +# | ENTRY POINT | +# '-------------' +if __name__ == '__main__': + # ,--------------, + # | Start timers | + # '--------------' + U.tic() + t1 = time.perf_counter() + + # ,------, + # | LOGO | + # '------' + U.print_logo() + + # ,-----------------, + # | Present options | + # '-----------------' + parser = argparse.ArgumentParser() + parser.add_argument('-v', '--version', help='Displays current package version.', action='store_true') + parser.add_argument('-mg', help="Given a path, return GRS metadata", default='local', type=str) + args = parser.parse_args().__dict__ # Converts the input arguments from Namespace() to dict + + print('User Input Arguments:') + for key in args: + print(f'{key}: {args[key]}') + + # ,----------------------------------------, + # | Automation pipelines class declaration | + # '----------------------------------------' + gpk_pipe = A.Pipelines() + + # ,---------------, + # | Treat options | + # '---------------' + if args['version']: + print(f'GET-Pak version: {getpak.__version__}') + elif args['mg']: + gpk_pipe.get_grs_metadict(args['mg']) + else: + print('No commands supplied, exiting.\n') + + # ,------------------------------, + # | End timers and report to log | + # '------------------------------' + t_hour, t_min, t_sec, _ = U.tac() + t2 = time.perf_counter() + final_msg_1 = f'Finished in {round(t2 - t1, 2)} second(s).' + final_msg_2 = f'Elapsed execution time: {t_hour}h : {t_min}m : {t_sec}s' + print(final_msg_1) + print(final_msg_2) + # ,-----, + # | END | + # '-----' diff --git a/requirements.txt b/requirements.txt index a199642..705602d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ +gdal fiona matplotlib netcdf4 @@ -5,4 +6,5 @@ numpy pandas rasterio rasterstats -scikit_learn \ No newline at end of file +scikit_learn +xarray \ No newline at end of file diff --git a/setup.py b/setup.py index 0303651..34149fe 100644 --- a/setup.py +++ b/setup.py @@ -1,38 +1,38 @@ -import getpak -from setuptools import setup, find_packages - -__version__ = getpak.__version__ -__package__ = 'getpak' -short_description = 'Sentinel-2 and 3 raster and vector manipulation and validation tools.' - -setup( - name=__package__, - version=__version__, - url='https://github.com/daviguima/get-pak', - packages=find_packages(), - classifiers=[ - "Programming Language :: Python :: 3", - "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", - "Operating System :: OS Independent", - ], - package_data={ - '': ['*.json'], - 'getpak': ['getpak/data/*'] - }, - include_package_data=True, - - license='GPLv3', - author='David Guimaraes', - author_email='dvdgmf@gmail.com', - description=short_description, - entry_points={ - 'console_scripts': ['getpak=main:main'], - }, - install_requires=[ - 'scikit_learn', - 'matplotlib', - 'numpy', - 'pandas', - 'rasterstats' - ] - ) +import getpak +from setuptools import setup, find_packages + +__version__ = getpak.__version__ +__package__ = 'getpak' +short_description = 'Sentinel-2 and 3 raster and vector manipulation and validation tools.' + +setup( + name=__package__, + version=__version__, + url='https://github.com/daviguima/get-pak', + packages=find_packages(), + classifiers=[ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Operating System :: OS Independent", + ], + package_data={ + '': ['*.json'], + 'getpak': ['getpak/data/*'] + }, + include_package_data=True, + + license='GPLv3', + author='David Guimaraes', + author_email='dvdgmf@gmail.com', + description=short_description, + entry_points={ + 'console_scripts': ['getpak=main:main'], + }, + install_requires=[ + 'scikit_learn', + 'matplotlib', + 'numpy', + 'pandas', + 'rasterstats' + ] + )