Skip to content
Snippets Groups Projects
Commit bb44a9ea authored by TGermain's avatar TGermain
Browse files

updated indices formulas for sen2cor baseline 4&5

parent 89f746c0
No related branches found
No related tags found
No related merge requests found
Pipeline #1769 passed
Stage: test
    Stage: deploy
      Hide whitespace changes
      Inline Side-by-side
      sen2chain/indices_functions.py
      + 69
      − 15
      View file @ bb44a9ea
      ......@@ -14,6 +14,7 @@ from rasterio.warp import reproject, Resampling
      import numpy as np
      from typing import Union
      from osgeo import gdal
      from packaging import version
      # import gdal
      ......@@ -93,6 +94,7 @@ logging.basicConfig(level=logging.INFO)
      def create_raw_ndr(
      baseline,
      b1_path: Union[str, pathlib.PosixPath],
      b2_path: Union[str, pathlib.PosixPath],
      out_path: Union[str, pathlib.PosixPath] = "./raw_ndr.tif",
      ......@@ -135,8 +137,14 @@ def create_raw_ndr(
      ),
      resampling=Resampling.bilinear,
      ).astype(np.float32)
      ndr = ((b1 - b2) / (b1 + b2) * 10000).astype(np.int16)
      ndr_masked = np.where(b1 != 0, ndr, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      ndr = (((b1-1000) - (b2-1000)) / ((b1-1000) + (b2-1000)) * 10000).astype(np.int16)
      ndr_masked = np.where(b1 != 0, ndr, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      ndr = ((b1 - b2) / (b1 + b2) * 10000).astype(np.int16)
      ndr_masked = np.where(b1 != 0, ndr, 32767)
      b1_profile.update(
      driver="Gtiff",
      ......@@ -296,6 +304,7 @@ def create_raw_ndr(
      def create_raw_ireci(
      baseline,
      b1_path: Union[str, pathlib.PosixPath],
      b2_path: Union[str, pathlib.PosixPath],
      b3_path: Union[str, pathlib.PosixPath],
      ......@@ -366,9 +375,14 @@ def create_raw_ireci(
      ),
      resampling=Resampling.bilinear,
      ).astype(np.float32)
      ireci = (b4 * (b1 - b2) / b3).astype(np.int16)
      ireci_masked = np.where(b1 != 0, ireci, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      ireci = ((b4-1000) * ((b1-1000) - (b2-1000)) / (b3-1000)).astype(np.int16)
      ireci_masked = np.where(b1 != 0, ireci, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      ireci = (b4 * (b1 - b2) / b3).astype(np.int16)
      ireci_masked = np.where(b1 != 0, ireci, 32767)
      b2_profile.update(
      driver="Gtiff",
      ......@@ -386,6 +400,7 @@ def create_raw_ireci(
      def create_raw_bigr(
      baseline,
      red_path: Union[str, pathlib.PosixPath],
      green_path: Union[str, pathlib.PosixPath],
      out_path: Union[str, pathlib.PosixPath] = "./raw_bigr.tif",
      ......@@ -408,8 +423,15 @@ def create_raw_bigr(
      np.seterr(
      divide="ignore", invalid="ignore"
      ) # ignore warnings when dividing by zero
      bigr = ((((green) ** 2 + (red) ** 2) / 2) ** 0.5).astype(np.int16)
      bigr_masked = np.where(red != 0, bigr, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      bigr = ((((green-1000) ** 2 + (red-1000) ** 2) / 2) ** 0.5).astype(np.int16)
      bigr_masked = np.where(red != 0, bigr, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      bigr = ((((green) ** 2 + (red) ** 2) / 2) ** 0.5).astype(np.int16)
      bigr_masked = np.where(red != 0, bigr, 32767)
      red_profile.update(
      driver="Gtiff",
      ......@@ -427,6 +449,7 @@ def create_raw_bigr(
      def create_raw_birnir(
      baseline,
      red_path: Union[str, pathlib.PosixPath],
      nir_path: Union[str, pathlib.PosixPath],
      out_path: Union[str, pathlib.PosixPath] = "./raw_birnir.tif",
      ......@@ -449,8 +472,15 @@ def create_raw_birnir(
      np.seterr(
      divide="ignore", invalid="ignore"
      ) # ignore warnings when dividing by zero
      birnir = ((((nir) ** 2 + (red) ** 2) / 2) ** 0.5).astype(np.int16)
      birnir_masked = np.where(red != 0, birnir, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      birnir = ((((nir-1000) ** 2 + (red-1000) ** 2) / 2) ** 0.5).astype(np.int16)
      birnir_masked = np.where(red != 0, birnir, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      birnir = ((((nir) ** 2 + (red) ** 2) / 2) ** 0.5).astype(np.int16)
      birnir_masked = np.where(red != 0, birnir, 32767)
      red_profile.update(
      driver="Gtiff",
      ......@@ -468,6 +498,7 @@ def create_raw_birnir(
      def create_raw_bibg(
      baseline,
      blue_path: Union[str, pathlib.PosixPath],
      green_path: Union[str, pathlib.PosixPath],
      out_path: Union[str, pathlib.PosixPath] = "./raw_bibg.tif",
      ......@@ -490,8 +521,15 @@ def create_raw_bibg(
      np.seterr(
      divide="ignore", invalid="ignore"
      ) # ignore warnings when dividing by zero
      bibg = ((((green) ** 2 + (blue) ** 2) / 2) ** 0.5).astype(np.int16)
      bibg_masked = np.where(blue != 0, bibg, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      bibg = ((((green-1000) ** 2 + (blue-1000) ** 2) / 2) ** 0.5).astype(np.int16)
      bibg_masked = np.where(blue != 0, bibg, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      bibg = ((((green) ** 2 + (blue) ** 2) / 2) ** 0.5).astype(np.int16)
      bibg_masked = np.where(blue != 0, bibg, 32767)
      blue_profile.update(
      driver="Gtiff",
      ......@@ -510,6 +548,7 @@ def create_raw_bibg(
      # A faire
      def create_raw_bi(
      baseline,
      b1_path: Union[str, pathlib.PosixPath],
      b2_path: Union[str, pathlib.PosixPath],
      out_path: Union[str, pathlib.PosixPath] = "./raw_bi.tif",
      ......@@ -532,8 +571,15 @@ def create_raw_bi(
      np.seterr(
      divide="ignore", invalid="ignore"
      ) # ignore warnings when dividing by zero
      bigr = ((((green) ** 2 + (red) ** 2) / 2) ** 0.5).astype(np.int16)
      bigr_masked = np.where(red != 0, bigr, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      bigr = ((((green-1000) ** 2 + (red-1000) ** 2) / 2) ** 0.5).astype(np.int16)
      bigr_masked = np.where(red != 0, bigr, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      bigr = ((((green) ** 2 + (red) ** 2) / 2) ** 0.5).astype(np.int16)
      bigr_masked = np.where(red != 0, bigr, 32767)
      red_profile.update(
      driver="Gtiff",
      ......@@ -550,6 +596,7 @@ def create_raw_bi(
      return Path(str(out_path)).absolute
      def create_raw_evi(
      baseline,
      blue_path: Union[str, pathlib.PosixPath],
      red_path: Union[str, pathlib.PosixPath],
      nir_path: Union[str, pathlib.PosixPath],
      ......@@ -573,8 +620,15 @@ def create_raw_evi(
      np.seterr(
      divide="ignore", invalid="ignore"
      ) # ignore warnings when dividing by zero
      evi = (2.5 * ((nir/10000 - red/10000)) / ((nir/10000 + 6.0 * red/10000- 7.5 * blue/10000)+ 1)*10000).astype(np.int32)
      evi_masked = np.where(red != 0, evi, 32767)
      if version.parse(baseline) >= version.parse("4.0"):
      logger.info('Baseline >= 4.0 detected, changing DN calculation')
      evi = (2.5 * (((nir-1000)/10000) - ((red-1000)/10000)) / (((nir-1000)/10000) + 6.0 * ((red-1000)/10000) - 7.5 * ((blue-1000)/10000) + 1)*10000).astype(np.int32)
      evi_masked = np.where(red != 0, evi, 32767)
      else:
      logger.info('Baseline <= 4.0 detected, sticking with original DN calculation')
      evi = (2.5 * ((nir/10000 - red/10000)) / ((nir/10000 + 6.0 * red/10000 - 7.5 * blue/10000) + 1)*10000).astype(np.int32)
      evi_masked = np.where(red != 0, evi, 32767)
      red_profile.update(
      driver="Gtiff",
      ......
      0% or .
      You are about to add 0 people to the discussion. Proceed with caution.
      Finish editing this message first!
      Please register or to comment