Commit 22f9fbd5 authored by Karl Segl's avatar Karl Segl
Browse files

Merge branch 'master' into 'feature/add_deadpixelcorrection'

# Conflicts:
#   enpt/execution/controller.py
#   enpt/model/images.py
#   requirements.txt
#   setup.py
#   tests/gitlab_CI_docker/context/environment_enpt.yml
parents 7b985872 f68cbe8f
Pipeline #3002 failed with stage
in 16 minutes and 50 seconds
......@@ -13,6 +13,16 @@ test_enpt:
- source /root/miniconda3/bin/activate enpt
- export GDAL_DATA=/root/miniconda3/envs/enpt/share/gdal
- export PYTHONPATH=$PYTHONPATH:/root # /root <- here are the sicor tables
# update sicor
# - conda install -y -q -c conda-forge basemap
# - rm -rf context/sicor
# - git clone https://gitext.gfz-potsdam.de/EnMAP/sicor.git ./context/sicor
# - cd ./context/sicor
# - make download-tables
# - python setup.py install
# - cd ../../
# run nosetests
- make nosetests # test are called here
# create the docs
......@@ -36,6 +46,7 @@ test_styles:
- make lint
artifacts:
paths:
- tests/data/test_outputs/*.log
- tests/linting/flake8.log
- tests/linting/pycodestyle.log
- tests/linting/pydocstyle.log
......@@ -48,16 +59,28 @@ test_enpt_install:
- source /root/miniconda3/bin/activate
- conda create -y -q --name enpt_test python=3
- source activate enpt_test
# install some dependencies that cause trouble when installed via pip
- conda install -y -c conda-forge scipy # scikit-image, matplotlib
# install not pip-installable deps of geoarray
- conda install -y -c conda-forge numpy scikit-image matplotlib pandas gdal rasterio pyproj basemap shapely
- conda install -y -c conda-forge 'icu=58.*' lxml # fixes bug for conda-forge gdal build
# install sicor
- conda install -y -q -c conda-forge pygrib h5py pytables pyfftw numba llvmlite
- rm -rf context/sicor
- git clone https://gitext.gfz-potsdam.de/EnMAP/sicor.git ./context/sicor
- cd ./context/sicor
- make install
- cd ../../
# install enpt
- make install
- cd ..
- pwd
- ls
# test importability
- python -c "import enpt; print(enpt)"
- python -c "from enpt.model.images import EnMAPL1Product_SensorGeo"
......
......@@ -6,7 +6,7 @@ EnPT - EnMAP Processing Tools
Please check the documentation_ for usage and in depth information.
Licence
License
-------
Free software: GNU General Public License v3
......@@ -25,7 +25,7 @@ See also the latest coverage_ report and the nosetests_ HTML report.
Credits
---------
-------
This software was funded from BMBF under EnMAP ...
......
......@@ -11,7 +11,7 @@ import warnings
from ..options.config import EnPTConfig
from ..io.reader import L1B_Reader
from ..processors.radiometric_transform import Radiometric_Transformer
from ..processors import Radiometric_Transformer
from ..model.images import EnMAPL1Product_SensorGeo
......@@ -57,7 +57,7 @@ class EnPT_Controller(object):
if not os.path.isdir(path_enmap_image) and \
not (os.path.exists(path_enmap_image) and path_enmap_image.endswith('.zip')):
raise ValueError("The parameter 'path_enmap_image' must be a directory or the path to an existing zip "
"archive.")
"archive. Received %s." % path_enmap_image)
# extract L1B image archive if needed
if path_enmap_image.endswith('.zip'):
......@@ -79,9 +79,12 @@ class EnPT_Controller(object):
# run transformation to TOARef
self.L1_obj = RT.transform_TOARad2TOARef(self.L1_obj)
def run_geometry_processor(self):
pass
def run_atmospheric_correction(self):
"""Run atmospheric correction only."""
pass
self.L1_obj.run_AC()
def write_output(self):
if self.cfg.output_dir:
......
......@@ -49,7 +49,7 @@ class _EnMAP_Image(object):
self.paths = SimpleNamespace()
# protected attributes
self._arr = None
self._data = None
self._mask_nodata = None
self._mask_clouds = None
self._mask_clouds_confidence = None
......@@ -63,7 +63,7 @@ class _EnMAP_Image(object):
self.basename = ''
@property
def data(self):
def data(self) -> GeoArray:
"""Return the actual EnMAP image data.
Bundled with all the corresponding metadata.
......@@ -107,26 +107,25 @@ class _EnMAP_Image(object):
:return: instance of geoarray.GeoArray
"""
# TODO this must return a subset if self.subset is not None
return self._arr
return self._data
@data.setter
def data(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
# TODO this must be able to handle subset inputs in tiled processing
if isinstance(geoArr_initArgs[0], np.ndarray):
self._arr = GeoArray(geoArr_initArgs[0], geotransform=self._arr.gt, projection=self._arr.prj)
self._data = GeoArray(geoArr_initArgs[0], geotransform=self._data.gt, projection=self._data.prj)
else:
self._arr = GeoArray(*geoArr_initArgs)
self._data = GeoArray(*geoArr_initArgs)
else:
del self.data
@data.deleter
def data(self):
self._arr = None
self._data = None
@property
def mask_clouds(self):
def mask_clouds(self) -> GeoArray:
"""Return the cloud mask.
Bundled with all the corresponding metadata.
......@@ -157,7 +156,7 @@ class _EnMAP_Image(object):
self._mask_clouds = None
@property
def mask_clouds_confidence(self):
def mask_clouds_confidence(self) -> GeoArray:
"""Return pixelwise information on the cloud mask confidence.
Bundled with all the corresponding metadata.
......@@ -192,7 +191,7 @@ class _EnMAP_Image(object):
self._mask_clouds_confidence = None
@property
def dem(self):
def dem(self) -> GeoArray:
"""Return an SRTM DEM in the exact dimension an pixel grid of self.arr.
:return: geoarray.GeoArray
......@@ -222,7 +221,7 @@ class _EnMAP_Image(object):
self._dem = None
@property
def ac_errors(self):
def ac_errors(self) -> GeoArray:
"""Return error information calculated by the atmospheric correction.
:return: geoarray.GeoArray
......@@ -253,7 +252,7 @@ class _EnMAP_Image(object):
self._ac_errors = None
@property
def deadpixelmap(self):
def deadpixelmap(self) -> GeoArray:
"""Return the dead pixel map.
Bundled with all the corresponding metadata. Dimensions: (bands x columns).
......@@ -517,6 +516,11 @@ class EnMAPL1Product_SensorGeo(object):
"""Correct dead pixels of both detectors."""
self.vnir.correct_dead_pixels()
self.swir.correct_dead_pixels()
def run_AC(self):
from ..processors import AtmosphericCorrector
AC = AtmosphericCorrector(config=self.cfg)
AC.run_ac(self)
def save(self, outdir: str, suffix="") -> str:
"""Save this product to disk using almost the same format as for reading.
......@@ -599,7 +603,7 @@ class EnMAP_Detector_MapGeo(_EnMAP_Image):
super(EnMAP_Detector_MapGeo, self).__init__()
@property
def mask_nodata(self):
def mask_nodata(self) -> GeoArray:
"""Return the no data mask.
Bundled with all the corresponding metadata.
......@@ -633,7 +637,7 @@ class EnMAP_Detector_MapGeo(_EnMAP_Image):
def mask_nodata(self):
self._mask_nodata = None
def calc_mask_nodata(self, fromBand=None, overwrite=False):
def calc_mask_nodata(self, fromBand=None, overwrite=False) -> GeoArray:
"""Calculate a no data mask with (values: 0=nodata; 1=data).
:param fromBand: <int> index of the band to be used (if None, all bands are used)
......
......@@ -65,6 +65,7 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
self.geom_view_azimuth = None # type: float # viewing azimuth angle
self.geom_sun_zenith = None # type: float # sun zenith angle
self.geom_sun_azimuth = None # type: float # sun azimuth angle
self.mu_sun = None # type: float # needed by SICOR for TOARad > TOARef conversion
self.lat_UL_UR_LL_LR = None # type: list # latitude coordinates for UL, UR, LL, LR
self.lon_UL_UR_LL_LR = None # type: list # longitude coordinates for UL, UR, LL, LR
self.lats = None # type: np.ndarray # 2D array of latitude coordinates according to given lon/lat sampling
......@@ -106,6 +107,7 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
xml.findall("%s/illumination_geometry/zenith_angle" % lbl)[0].text.split()[0])
self.geom_sun_azimuth = np.float(
xml.findall("%s/illumination_geometry/azimuth_angle" % lbl)[0].text.split()[0])
self.mu_sun = np.cos(np.deg2rad(self.geom_sun_zenith))
self.lat_UL_UR_LL_LR = \
[float(xml.findall("%s/geometry/bounding_box/%s_northing" % (lbl, corner))[0].text.split()[0])
for corner in ("UL", "UR", "LL", "LR")]
......
......@@ -17,6 +17,8 @@ from collections import OrderedDict, Mapping
import numpy as np
from multiprocessing import cpu_count
import sicor
from .options_schema import \
enpt_schema_input, \
enpt_schema_config_output, \
......@@ -31,6 +33,14 @@ path_enptlib = os.path.dirname(pkgutil.get_loader("enpt").path)
path_options_default = os.path.join(path_enptlib, 'options', 'options_default.json')
config_for_testing = dict(
path_l1b_enmap_image=os.path.abspath(
os.path.join(path_enptlib, '..', 'tests', 'data', 'EnMAP_Level_1B', 'AlpineTest1_CWV2_SM0.zip')),
log_level='DEBUG',
output_dir=os.path.join(path_enptlib, '..', 'tests', 'data', 'test_outputs')
)
class EnPTConfig(object):
def __init__(self, json_config='', **user_opts):
"""Create a job configuration.
......@@ -72,7 +82,7 @@ class EnPTConfig(object):
# output options #
##################
self.output_dir = self.absPath(gp('output_dir'))
self.output_dir = self.absPath(gp('output_dir', fallback=os.path.abspath(os.path.curdir)))
###########################
# processor configuration #
......@@ -81,6 +91,7 @@ class EnPTConfig(object):
# toa_ref
self.path_earthSunDist = self.absPath(gp('path_earthSunDist'))
self.path_solar_irr = self.absPath(gp('path_solar_irr'))
self.scale_factor_toa_ref = gp('scale_factor_toa_ref')
# geometry
self.enable_keystone_correction = gp('enable_keystone_correction')
......@@ -88,8 +99,10 @@ class EnPTConfig(object):
self.path_reference_image = gp('path_reference_image')
# atmospheric_correction
self.sicor_cache_dir = gp('sicor_cache_dir', fallback=sicor.__path__[0])
self.auto_download_ecmwf = gp('auto_download_ecmwf')
self.enable_cloud_screening = gp('enable_cloud_screening')
self.scale_factor_boa_ref = gp('scale_factor_boa_ref'),
# smile
self.run_smile_P = gp('run_smile_P')
......
......@@ -16,7 +16,8 @@
"processors": {
"toa_ref": {
"path_earthSunDist": "./resources/earth_sun_distance/Earth_Sun_distances_per_day_edited__1980_2030.csv", /*input path of the earth sun distance model*/
"path_solar_irr": "./resources/solar_irradiance/SUNp1fontenla__350-2500nm_@0.1nm_converted.txt" /*input path of the solar irradiance model*/
"path_solar_irr": "./resources/solar_irradiance/SUNp1fontenla__350-2500nm_@0.1nm_converted.txt", /*input path of the solar irradiance model*/
"scale_factor_toa_ref": 10000 /*scale factor to be applied to TOA reflectance result*/
},
"geometry": {
......@@ -26,8 +27,12 @@
},
"atmospheric_correction": {
"sicor_cache_dir": "", /*directory to be used to stored sicor cache files
NOTE: SICOR stores intermediate results there that need to computed only once
for atmospheric correction of multiple EnMAP images. (default: 'auto')*/
"auto_download_ecmwf": false,
"enable_cloud_screening": false
"enable_cloud_screening": false,
"scale_factor_boa_ref": 10000 /*scale factor to be applied to BOA reflectance result*/
},
"smile": {
......
......@@ -41,6 +41,7 @@ enpt_schema_input = dict(
atmospheric_correction=dict(
type='dict', required=False,
schema=dict(
sicor_cache_dir=dict(type='string', required=False),
auto_download_ecmwf=dict(type='boolean', required=False),
enable_cloud_screening=dict(type='boolean', required=False),
)),
......@@ -85,6 +86,7 @@ parameter_mapping = dict(
# processors > toa_ref
path_earthSunDist=('processors', 'toa_ref', 'path_earthSunDist'),
path_solar_irr=('processors', 'toa_ref', 'path_solar_irr'),
scale_factor_toa_ref=('processors', 'toa_ref', 'scale_factor_toa_ref'),
# processors > geometry
enable_keystone_correction=('processors', 'geometry', 'enable_keystone_correction'),
......@@ -92,8 +94,10 @@ parameter_mapping = dict(
path_reference_image=('processors', 'geometry', 'path_reference_image'),
# processors > atmospheric_correction
sicor_cache_dir=('processors', 'atmospheric_correction', 'sicor_cache_dir'),
auto_download_ecmwf=('processors', 'atmospheric_correction', 'auto_download_ecmwf'),
enable_cloud_screening=('processors', 'atmospheric_correction', 'enable_cloud_screening'),
scale_factor_boa_ref=('processors', 'atmospheric_correction', 'scale_factor_boa_ref'),
# processors > smile
run_smile_P=('processors', 'smile', 'run_processor'),
......
# -*- coding: utf-8 -*-
"""EnPT 'processors' module containing all EnPT processor sub-modules."""
from .radiometric_transform.radiometric_transform import Radiometric_Transformer
from .atmospheric_correction.atmospheric_correction import AtmosphericCorrector
__all__ = [
"Radiometric_Transformer",
"AtmosphericCorrector"
]
# -*- coding: utf-8 -*-
"""EnPT 'atmospheric correction module."""
from .atmospheric_correction import AtmosphericCorrector
__all__ = ['AtmosphericCorrector']
......@@ -3,3 +3,69 @@
Performs the atmospheric correction of EnMAP L1B data.
"""
import pprint
import numpy as np
from sicor.sicor_enmap import sicor_ac_enmap
from sicor.options import get_options as get_ac_options
from ...model.images import EnMAPL1Product_SensorGeo
from ...options.config import EnPTConfig
from ...utils.path_generator import get_path_ac_options
class AtmosphericCorrector(object):
"""Class for performing atmospheric correction of EnMAP L1 images using SICOR."""
def __init__(self, config: EnPTConfig=None):
"""Create an instance of AtmosphericCorrector."""
self.cfg = config
@staticmethod
def get_ac_options(buffer_dir):
path_opts = get_path_ac_options()
try:
options = get_ac_options(path_opts)
# adjust options
options['EnMAP']['buffer_dir'] = buffer_dir
for vv in options["RTFO"].values():
vv["hash_formats"] = dict(spr='%.0f',
coz='%.0f,',
cwv='%.0f,',
tmp='%0f,',
tau_a='%.2f,',
vza='%.0f,')
options["ECMWF"]["path_db"] = "./ecmwf"
return options
except FileNotFoundError:
raise FileNotFoundError('Could not locate options file for atmospheric correction at %s.' % path_opts)
def run_ac(self, enmap_ImageL1: EnMAPL1Product_SensorGeo) -> EnMAPL1Product_SensorGeo:
options = self.get_ac_options(buffer_dir=self.cfg.sicor_cache_dir)
enmap_ImageL1.logger.debug('AC options: \n' + pprint.pformat(options))
# run AC
enmap_ImageL1.logger.info("Starting atmospheric correction for VNIR and SWIR detector. "
"Source radiometric unit code is '%s'." % enmap_ImageL1.meta.vnir.unitcode)
enmap_l2a_sens_geo, state, cwv_map, ch4_map = sicor_ac_enmap(enmap_l1b=enmap_ImageL1, options=options,
logger=enmap_ImageL1.logger, debug=True)
# join results
enmap_ImageL1.logger.info('Joining results of atmospheric correction.')
for in_detector, out_detector in zip([enmap_ImageL1.vnir, enmap_ImageL1.swir],
[enmap_l2a_sens_geo.vnir, enmap_l2a_sens_geo.swir]):
in_detector.data = (out_detector.data[:] * self.cfg.scale_factor_boa_ref).astype(np.int16)
# NOTE: geotransform and projection are missing due to sensor geometry
del in_detector.data_l2a # FIXME sicor sets data_l2a to float array -> not needed
del in_detector.unit # FIXME sicor sets unit to '1' -> not needed
in_detector.detector_meta.unit = '0-%d' % self.cfg.scale_factor_boa_ref
in_detector.detector_meta.unitcode = 'BOARef'
return enmap_ImageL1
# -*- coding: utf-8 -*-
"""EnPT 'radiometric transform' module containing eveything related to radiometric transformations."""
from .radiometric_transform import Radiometric_Transformer
__all__ = ['Radiometric_Transformer']
......@@ -21,8 +21,7 @@ class Radiometric_Transformer(object):
self.solarIrr = config.path_solar_irr # path of model for solar irradiance
self.earthSunDist = config.path_earthSunDist # path of model for earth sun distance
@staticmethod
def transform_TOARad2TOARef(enmap_ImageL1: EnMAPL1Product_SensorGeo, scale_factor: int=10000):
def transform_TOARad2TOARef(self, enmap_ImageL1: EnMAPL1Product_SensorGeo):
"""Transform top-of-atmosphere radiance to top-of-atmosphere reflectance.
NOTE: The following formula is used:
......@@ -30,7 +29,6 @@ class Radiometric_Transformer(object):
(solIrr * math.cos(zenithAngleDeg))
:param enmap_ImageL1: instance of the class 'EnMAPL1Product_ImGeo'
:param scale_factor: scale factor to be applied to TOA reflectance result
:return:
"""
for detectorName in enmap_ImageL1.detector_attrNames:
......@@ -41,7 +39,7 @@ class Radiometric_Transformer(object):
# compute TOA reflectance
constant = \
scale_factor * math.pi * enmap_ImageL1.meta.earthSunDist ** 2 / \
self.cfg.scale_factor_toa_ref * math.pi * enmap_ImageL1.meta.earthSunDist ** 2 / \
(math.cos(math.radians(detector.detector_meta.geom_sun_zenith)))
solIrr = np.array([detector.detector_meta.solar_irrad[band] for band in detector.detector_meta.srf.bands])\
.reshape(1, 1, detector.data.bands)
......@@ -49,7 +47,7 @@ class Radiometric_Transformer(object):
# update EnMAP image
detector.data = toaRef
detector.detector_meta.unit = '0-%d' % scale_factor
detector.detector_meta.unit = '0-%d' % self.cfg.scale_factor_toa_ref
detector.detector_meta.unitcode = 'TOARef'
return enmap_ImageL1
......@@ -48,3 +48,11 @@ class PathGenL1BProduct(object):
def _find_in_metaxml(self, expression):
return self.xml.findall(expression)[0].text.replace("\n", "").strip()
def get_path_ac_options() -> str:
"""Returns the path of the options json file needed for atmospheric correction."""
from sicor import options
path_ac = os.path.join(os.path.dirname(options.__file__), 'sicor_enmap_user_options.json')
return path_ac
......@@ -6,3 +6,4 @@ cerberus
jsmin
matplotlib
tqdm
git+https://gitext.gfz-potsdam.de/EnMAP/sicor.git
\ No newline at end of file
......@@ -6,5 +6,3 @@ flake8==2.6.0
tox==2.3.1
coverage==4.1
Sphinx==1.4.8
......@@ -15,6 +15,7 @@ with open("enpt/version.py") as version_file:
requirements = [ # put package requirements here
'numpy', 'scipy', 'geoarray>=0.7.15', 'spectral>=0.16', 'cerberus', 'jsmin', 'matplotlib', 'tqdm'
# 'sicor', # pip install git+https://gitext.gfz-potsdam.de/EnMAP/sicor.git
]
test_requirements = ['coverage', 'nose', 'nose-htmloutput', 'rednose']
......
# -*- coding: utf-8 -*-
import os
enptRepo_rootpath = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
config_for_testing = dict(
path_l1b_enmap_image=os.path.join(enptRepo_rootpath, 'tests', 'data', 'EnMAP_Level_1B', 'AlpineTest1_CWV2_SM0.zip'),
output_dir=os.path.join(enptRepo_rootpath, 'tests', 'data', 'test_outputs')
)
......@@ -2,9 +2,19 @@
context_dir="./context"
dockerfile="enpt_ci.docker"
tag="enpt_ci:latest"
tag="enpt_ci:0.4.0b4"
gitlab_runner="enpt_gitlab_CI_runner"
# get sicor project
rm -rf context/sicor
# git clone https://gitext.gfz-potsdam.de/EnMAP/sicor.git ./context/sicor
git clone https://gitext.gfz-potsdam.de/EnMAP/sicor.git --branch feature/improve_enmap --single-branch ./context/sicor
# download sicor cache (fastens SICOR CI tests a lot, but cache needs to be updated manually using a local sicor repo:
# 1. clone a fresh copy of sicor or delete sicor/sicor/aerosol_0_ch4_34d3778719cc87188787de09bb8f870d16050078.pkl.zip
# 2. run a sicor test including sicor_ac or enmap_ac (recreates cache file) -> upload newly created cache file
# wget http://ouo.io/uCQxof -P ./context/
echo "#### Build runner docker image"
sudo docker rmi ${tag}
sudo docker build -f ${context_dir}/${dockerfile} -m 20G -t ${tag} ${context_dir}
......
......@@ -4,12 +4,12 @@ FROM centos:7
RUN yum update -y && \
yum install -y wget vim bzip2 gcc gcc-c++ make libgl1-mesa-glx mesa-libGL qt5-qtbase-gui git nano tree gdb texlive
# change version numbers for future upgrades
# ENV miniconda_dl 'Miniconda3-latest-Linux-x86_64.sh'
# change version numbers for future upgrades # currently Anaconda 4.5.1 (worked fine)
ENV miniconda_dl 'Miniconda3-latest-Linux-x86_64.sh'
# use miniconda 4.3.31 as workaround for "IsADirectoryError(21, 'Is a directory')" with version 4.4.10 (currently latest)
ENV miniconda_dl 'Miniconda3-4.3.31-Linux-x86_64.sh'
# ENV miniconda_dl 'Miniconda3-4.3.31-Linux-x86_64.sh'
ENV envconfig 'environment_enpt.yml'
ENV git_lfs_v='2.1.1'
ENV git_lfs_v='2.4.1'
RUN /bin/bash -i -c "cd /root; wget https://repo.continuum.io/miniconda/$miniconda_dl ; \
bash -i /root/$miniconda_dl -b ; \
......@@ -30,3 +30,19 @@ RUN /bin/bash -i -c "wget https://github.com/git-lfs/git-lfs/releases/download/v
tar -zxvf git-lfs-linux-amd64-${git_lfs_v}.tar.gz; \
cd git-lfs-${git_lfs_v}; \
bash ./install.sh"
# copy sicor code to /tmp
COPY sicor /tmp/sicor
# install sicor (in pip development mode so that its root directory in /tmp/sicor matching with the subsequent COPY command)
RUN bash -i -c "source /root/miniconda3/bin/activate enpt; \
cd /tmp/sicor/ ; \
make clean ; \
make requirements ; \
make download-tables ; \
pip install -e . --no-cache-dir"
# copy sicor cache files to sicor root directory (speeds up SICOR CI tests because table subsets dont have to be created each time)
# -> sicor root directory is the default directory of these cache files if sicor_cache_dir is not set in EnPT options
COPY *.zip /tmp/sicor/sicor
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment