diff --git a/Docker/Dockerfile b/Docker/Dockerfile index 0d250c66..1367665d 100644 --- a/Docker/Dockerfile +++ b/Docker/Dockerfile @@ -1,8 +1,8 @@ FROM ubuntu:22.04 LABEL author="OPERA ADT" \ - description="RTC interface release" \ - version="interface_0.1" + description="RTC beta release" \ + version="0.2-beta" RUN apt-get -y update &&\ apt-get -y install curl zip make &&\ @@ -17,7 +17,7 @@ USER rtc_user ENV CONDA_PREFIX=/home/rtc_user/miniconda3 -#Install conda and +# Install conda and RTC WORKDIR /home/rtc_user RUN curl -sSL https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -o miniconda.sh &&\ bash miniconda.sh -b -p ${CONDA_PREFIX} &&\ @@ -25,60 +25,25 @@ RUN curl -sSL https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64 rm ${HOME}/miniconda.sh ENV PATH=${CONDA_PREFIX}/bin:${PATH} -ENV POSTFIX_ISCE3="rtc" -ENV ISCEHOME=/home/rtc_user - RUN ${CONDA_PREFIX}/bin/conda init bash -#Install ISCE3 from source code -RUN echo "Installing Conda environment for ISCE3" &&\ - mkdir -p ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/src &&\ - mkdir -p ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/build &&\ - mkdir -p ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/install &&\ - cd ${ISCEHOME} &&\ - curl -sSL https://github.com/gshiroma/isce3/archive/refs/heads/opera_rtc.zip -o isce3_opera_rtc.zip &&\ - unzip -d ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3} isce3_opera_rtc.zip &&\ - cd ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/src &&\ - ln -s ../isce3-opera_rtc isce3 &&\ - rm ${ISCEHOME}/isce3_opera_rtc.zip &&\ - cd ${HOME} &&\ - conda create --name isce3_rtc --file ./OPERA/RTC/Docker/specfile.txt &&\ - chmod -R 755 ${CONDA_PREFIX}/envs/isce3_rtc &&\ - cd ${ISCEHOME}/tools/isce_$POSTFIX_ISCE3/src - -SHELL ["conda", "run", "-n", "isce3_rtc", "/bin/bash", "-c"] +RUN conda create --name "RTC" --file /home/rtc_user/OPERA/RTC/Docker/specfile.txt +SHELL ["conda", "run", "-n", "RTC", "/bin/bash", "-c"] -ENV CC=${CONDA_PREFIX}/envs/isce3_rtc/bin/x86_64-conda-linux-gnu-gcc -ENV CXX=${CONDA_PREFIX}/envs/isce3_rtc/bin/x86_64-conda-linux-gnu-g++ -ENV PATH=$PATH:${ISCEHOME}/tools/isce_$POSTFIX_ISCE3/install/bin -ENV PYTHONPATH=$PYTHONPATH:${ISCEHOME}/tools/isce_$POSTFIX_ISCE3/install/packages -ENV LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/install/lib - -RUN cd ${ISCEHOME}/tools/isce_$POSTFIX_ISCE3/build &&\ - cmake -DCMAKE_INSTALL_PREFIX=${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/install -DWITH_CUDA=OFF ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/src/isce3 &&\ - make -j32 VERBOSE=ON &&\ - make install &&\ - rm -r ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/build &&\ - rm -r ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/isce3-opera_rtc &&\ - chmod -R 755 ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/install &&\ - cd /home/rtc_user/miniconda3/envs/isce3_rtc/lib &&\ - ln -s ${ISCEHOME}/tools/isce_${POSTFIX_ISCE3}/install/lib/libisce3.so &&\ - echo 'ISCE3 installation successful!' RUN echo "Installing OPERA s1-reader and RTC" &&\ - mkdir -p $HOME/OPERA &&\ cd $HOME/OPERA &&\ - curl -sSL https://github.com/seongsujeong/s1-reader/archive/refs/heads/correction_and_calibration.zip -o s1_reader.zip &&\ + curl -sSL https://github.com/opera-adt/s1-reader/archive/refs/tags/v0.1.3.zip -o s1_reader.zip &&\ unzip s1_reader.zip &&\ - chmod -R 755 s1-reader-correction_and_calibration &&\ - ln -s s1-reader-correction_and_calibration s1-reader &&\ + chmod -R 755 s1-reader-0.1.3 &&\ + ln -s s1-reader-0.1.3 s1-reader &&\ rm s1_reader.zip &&\ python -m pip install ./s1-reader &&\ - cd $HOME/OPERA &&\ - pip install ./RTC &&\ + python -m pip install ./RTC &&\ mkdir /home/rtc_user/scratch &&\ - echo 'conda activate isce3_rtc' >>/home/rtc_user/.bashrc + chmod -R 755 /home/rtc_user/scratch &&\ + echo 'conda activate RTC' >>/home/rtc_user/.bashrc WORKDIR /home/rtc_user/scratch -ENTRYPOINT ["conda", "run", "-n", "isce3_rtc","rtc_s1.py"] +ENTRYPOINT ["conda", "run", "--no-capture-output", "-n", "RTC","rtc_s1.py"] diff --git a/Docker/environment.yml b/Docker/environment.yml index a6ea2703..631722d8 100644 --- a/Docker/environment.yml +++ b/Docker/environment.yml @@ -9,9 +9,6 @@ dependencies: - gdal - gmock - gtest - - gcc_linux-64 < 10 - - gxx_linux-64 < 10 - - gfortran_linux-64 - h5py - hdf5 - libgcc-ng @@ -26,3 +23,4 @@ dependencies: - shapely - yamale - backoff + - isce3=0.9.0 diff --git a/Docker/requirements.txt b/Docker/requirements.txt index d7c3bad1..4b5fe591 100644 --- a/Docker/requirements.txt +++ b/Docker/requirements.txt @@ -5,3 +5,8 @@ ruamel.yaml scipy pytest requests +sqlite3 +pyproj +shapely +yamale + diff --git a/Docker/specfile.txt b/Docker/specfile.txt index ad1f6b94..3c516579 100644 --- a/Docker/specfile.txt +++ b/Docker/specfile.txt @@ -3,63 +3,55 @@ # platform: linux-64 @EXPLICIT https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/ca-certificates-2022.6.15-ha878542_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/ca-certificates-2022.12.7-ha878542_0.conda https://conda.anaconda.org/conda-forge/noarch/font-ttf-dejavu-sans-mono-2.37-hab24e00_0.tar.bz2 https://conda.anaconda.org/conda-forge/noarch/font-ttf-inconsolata-3.000-h77eed37_0.tar.bz2 https://conda.anaconda.org/conda-forge/noarch/font-ttf-source-code-pro-2.038-h77eed37_0.tar.bz2 https://conda.anaconda.org/conda-forge/noarch/font-ttf-ubuntu-0.83-hab24e00_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/kernel-headers_linux-64-2.6.32-he073ed8_15.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.36.1-hea4e1c9_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libgcc-devel_linux-64-9.5.0-h367e8d2_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libgfortran5-12.1.0-hdcd56e2_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-devel_linux-64-9.5.0-h367e8d2_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-12.1.0-ha89aaad_16.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.39-hcc3a1bd_1.conda +https://conda.anaconda.org/conda-forge/linux-64/libgfortran5-12.2.0-h337968e_19.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-12.2.0-h46fd767_19.tar.bz2 https://conda.anaconda.org/conda-forge/noarch/poppler-data-0.4.11-hd8ed1ab_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/tzdata-2022c-h191b570_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/python_abi-3.11-3_cp311.conda +https://conda.anaconda.org/conda-forge/noarch/tzdata-2022g-h191b570_0.conda https://conda.anaconda.org/conda-forge/noarch/fonts-conda-forge-1-0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libgfortran-ng-12.1.0-h69a702a_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libgomp-12.1.0-h8d9b700_16.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/sysroot_linux-64-2.12-he073ed8_15.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libgfortran-ng-12.2.0-h69a702a_19.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libgomp-12.2.0-h65d4601_19.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-2_gnu.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/binutils_impl_linux-64-2.36.1-h193b22a_2.tar.bz2 https://conda.anaconda.org/conda-forge/noarch/fonts-conda-ecosystem-1-0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/binutils_linux-64-2.36-hf3e587d_10.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-12.1.0-h8d9b700_16.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-12.2.0-h65d4601_19.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-h7f98852_4.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.18.1-h7f98852_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/eigen-3.4.0-h4bd325d_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/expat-2.4.8-h27087fc_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/fftw-3.3.10-nompi_ha7695d1_103.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/freexl-1.0.6-h7f98852_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/geos-3.11.0-h27087fc_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/expat-2.5.0-h27087fc_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/fftw-3.3.10-nompi_hf0379b8_106.conda +https://conda.anaconda.org/conda-forge/linux-64/freexl-1.0.6-h166bdaf_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/geos-3.11.1-h27087fc_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/gettext-0.21.1-h27087fc_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/giflib-5.2.1-h36c2ea0_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gtest-1.11.0-h924138e_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/icu-70.1-h27087fc_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/jpeg-9e-h166bdaf_2.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/json-c-0.16-hc379101_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/keyutils-1.6.1-h166bdaf_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/lerc-4.0.0-h27087fc_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libdeflate-1.13-h166bdaf_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libdeflate-1.14-h166bdaf_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20191231-he28a2e2_2.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libev-4.33-h516909a_1.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libffi-3.4.2-h7f98852_5.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libiconv-1.16-h516909a_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libiconv-1.17-h166bdaf_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libnsl-2.0.0-h7f98852_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libopenblas-0.3.21-pthreads_h78a6416_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libsanitizer-9.5.0-hf86b28c_16.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libopenblas-0.3.21-pthreads_h78a6416_3.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.32.1-h7f98852_1000.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libuv-1.44.2-h166bdaf_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libwebp-base-1.2.4-h166bdaf_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.2.12-h166bdaf_2.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.2.13-h166bdaf_4.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/lz4-c-1.9.3-h9c3ff4c_1.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.3-h27087fc_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/nspr-4.32-h9c3ff4c_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/openssl-3.0.5-h166bdaf_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pcre-8.45-h9c3ff4c_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda +https://conda.anaconda.org/conda-forge/linux-64/openssl-3.0.7-h0b41bf4_1.conda https://conda.anaconda.org/conda-forge/linux-64/pixman-0.40.0-h36c2ea0_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/pthread-stubs-0.4-h36c2ea0_1001.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/rhash-1.4.3-h166bdaf_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/snappy-1.1.9-hbd366e4_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/tzcode-2022c-h166bdaf_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/qhull-2020.2-h4bd325d_2.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/snappy-1.1.9-hbd366e4_2.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/tzcode-2022g-h166bdaf_0.conda https://conda.anaconda.org/conda-forge/linux-64/xorg-kbproto-1.0.7-h7f98852_1002.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xorg-libice-1.0.10-h7f98852_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xorg-libxau-1.0.9-h7f98852_0.tar.bz2 @@ -69,91 +61,72 @@ https://conda.anaconda.org/conda-forge/linux-64/xorg-xextproto-7.3.0-h7f98852_10 https://conda.anaconda.org/conda-forge/linux-64/xorg-xproto-7.0.31-h7f98852_1007.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xz-5.2.6-h166bdaf_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/yaml-0.2.5-h7f98852_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gcc_impl_linux-64-9.5.0-h6c5bc03_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gettext-0.19.8.1-h73d1719_1008.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gmock-1.11.0-h924138e_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/hdf4-4.2.15-h9772cbc_4.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/krb5-1.19.3-h08a2579_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libblas-3.9.0-16_linux64_openblas.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20191231-he28a2e2_2.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.47.0-hff17c54_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libpng-1.6.37-h753d276_4.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/librttopo-1.1.0-hf730bdb_11.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.39.2-h753d276_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libpng-1.6.39-h753d276_0.conda +https://conda.anaconda.org/conda-forge/linux-64/librttopo-1.1.0-ha49c73b_12.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.40.0-h753d276_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libssh2-1.10.0-hf14f497_3.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libxcb-1.13-h7f98852_1004.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.9.14-h22db469_4.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.10.3-h7463322_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libzip-1.9.2-hc929e4a_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/pcre2-10.40-hc3806b6_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/readline-8.1.2-h0f457ee_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/tk-8.6.12-h27826a3_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/xerces-c-3.2.3-h55805fa_5.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xorg-libsm-1.2.3-hd9c2040_1000.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/zlib-1.2.12-h166bdaf_2.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/zlib-1.2.13-h166bdaf_4.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.2-h6239696_4.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/blosc-1.21.1-h83bc5f7_3.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/boost-cpp-1.74.0-h75c5d50_8.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/freetype-2.12.1-hca18f0e_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gcc_linux-64-9.5.0-h4258300_10.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gfortran_impl_linux-64-9.5.0-h3c9b8b6_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gxx_impl_linux-64-9.5.0-h6c5bc03_16.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/krb5-1.19.3-h08a2579_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/blosc-1.21.2-hafa529b_0.conda +https://conda.anaconda.org/conda-forge/linux-64/boost-cpp-1.78.0-h75c5d50_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/freetype-2.12.1-hca18f0e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/hdf4-4.2.15-h9772cbc_5.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/libcblas-3.9.0-16_linux64_openblas.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libglib-2.72.1-h2d90d5f_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libcurl-7.86.0-h2283fc2_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libglib-2.74.1-h606061b_1.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/liblapack-3.9.0-16_linux64_openblas.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libtiff-4.4.0-h0e0dad5_3.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/sqlite-3.39.2-h4ff8645_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libpq-15.1-h67c24c5_1.conda +https://conda.anaconda.org/conda-forge/linux-64/libtiff-4.4.0-h55922b4_4.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libxslt-1.1.37-h873f0b0_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/nss-3.82-he02c5a1_0.conda +https://conda.anaconda.org/conda-forge/linux-64/python-3.11.0-ha86cf86_0_cpython.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/sqlite-3.40.0-h4ff8645_0.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xorg-libx11-1.7.2-h7f98852_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/fontconfig-2.14.0-h8e229c2_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gfortran_linux-64-9.5.0-hdb51d14_10.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gxx_linux-64-9.5.0-h43f449f_10.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/lcms2-2.12-hddcbb42_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libcurl-7.83.1-h2283fc2_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libkml-1.3.0-h238a007_1014.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libpq-14.5-he2d8382_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/nss-3.78-h2350873_0.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/backoff-2.2.1-pyhd8ed1ab_0.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/cached_property-1.5.2-pyha770c72_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/cfitsio-4.2.0-hd9d235c_0.conda +https://conda.anaconda.org/conda-forge/linux-64/curl-7.86.0-h2283fc2_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/fontconfig-2.14.1-hc2a2eb6_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/hdf5-1.12.2-nompi_h4df4325_100.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/lcms2-2.14-h6ed2654_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libkml-1.3.0-h37653c0_1015.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/lxml-4.9.1-py311hc4dbab1_1.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/numpy-1.23.5-py311h7d28db0_0.conda https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.0-h7d73246_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/python-3.9.13-h2660328_0_cpython.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/postgresql-15.1-ha105346_1.conda +https://conda.anaconda.org/conda-forge/linux-64/proj-9.1.0-h93bde94_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/pyyaml-6.0-py311hd4cff14_5.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/ruamel.yaml.clib-0.2.7-py311h2582759_1.conda +https://conda.anaconda.org/conda-forge/noarch/setuptools-65.5.1-pyhd8ed1ab_0.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/wheel-0.38.4-pyhd8ed1ab_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/xerces-c-3.2.4-h55805fa_1.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xorg-libxext-1.3.4-h7f98852_1.tar.bz2 https://conda.anaconda.org/conda-forge/linux-64/xorg-libxrender-0.9.10-h7f98852_1003.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/attrs-22.1.0-pyh71513ae_1.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/backoff-2.1.2-pyhd8ed1ab_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/cached_property-1.5.2-pyha770c72_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/cairo-1.16.0-ha61ee94_1012.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/cfitsio-4.1.0-hd9d235c_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/cmake-3.24.1-h5432695_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/curl-7.83.1-h2283fc2_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/hdf5-1.12.2-nompi_h4df4325_100.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/iniconfig-1.1.1-pyh9f0ad1d_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/postgresql-14.5-ha7cec9f_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/proj-9.0.1-h93bde94_1.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/py-1.11.0-pyh6c4a22f_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/pyparsing-3.0.9-pyhd8ed1ab_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/python_abi-3.9-2_cp39.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/setuptools-65.3.0-pyhd8ed1ab_1.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/tomli-2.0.1-pyhd8ed1ab_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/wheel-0.37.1-pyhd8ed1ab_0.tar.bz2 https://conda.anaconda.org/conda-forge/noarch/cached-property-1.5.2-hd8ed1ab_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/cython-0.29.32-py39h5a03fae_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/geotiff-1.7.1-h4fc65e6_3.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/kealib-1.4.15-ha7026e8_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libdap4-3.20.6-hd7c4107_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libnetcdf-4.8.1-nompi_h21705cb_104.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libspatialite-5.0.1-h38b5f51_18.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/numpy-1.23.2-py39hba7629e_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/packaging-21.3-pyhd8ed1ab_0.tar.bz2 -https://conda.anaconda.org/conda-forge/noarch/pip-22.2.2-pyhd8ed1ab_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pluggy-1.0.0-py39hf3d152e_3.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/poppler-22.04.0-h8b295ee_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pybind11-global-2.10.0-py39hf939315_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pyre-1.11.2-py39h8b9aba8_2.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pyyaml-6.0-py39hb9d737c_4.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/ruamel.yaml.clib-0.2.6-py39hb9d737c_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/tiledb-2.11.1-h3f4058f_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/h5py-3.7.0-nompi_py39hd51670d_101.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/libgdal-3.5.1-h84f2e82_5.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pybind11-2.10.0-py39hf939315_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/pytest-7.1.2-py39hf3d152e_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/ruamel.yaml-0.17.21-py39hb9d737c_1.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/scipy-1.9.0-py39h8ba3f38_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/shapely-1.8.4-py39h68ae834_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/cairo-1.16.0-ha61ee94_1014.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/geotiff-1.7.1-ha76d385_4.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/kealib-1.5.0-ha7026e8_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libnetcdf-4.8.1-nompi_h261ec11_106.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/libspatialite-5.0.1-h7c8129e_22.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/pip-22.3.1-pyhd8ed1ab_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/pyre-1.11.2-py311hbbb8f27_3.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/ruamel.yaml-0.17.21-py311hd4cff14_2.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/scipy-1.9.3-py311h69910c8_2.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/shapely-2.0.0-py311h0f577a2_0.conda +https://conda.anaconda.org/conda-forge/linux-64/tiledb-2.13.0-h3f4058f_0.conda https://conda.anaconda.org/conda-forge/noarch/yamale-4.0.4-pyh6c4a22f_0.tar.bz2 -https://conda.anaconda.org/conda-forge/linux-64/gdal-3.5.1-py39h5e8c895_5.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/h5py-3.7.0-nompi_py311hbe7f6d8_102.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/poppler-22.12.0-h92391eb_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libgdal-3.6.0-h4fca5ce_14.conda +https://conda.anaconda.org/conda-forge/linux-64/gdal-3.6.0-py311hadb6153_14.conda +https://conda.anaconda.org/conda-forge/linux-64/isce3-0.9.0-py311h0802494_0.conda \ No newline at end of file diff --git a/app/rtc_compare.py b/app/rtc_compare.py index cf362b87..d8994ea7 100644 --- a/app/rtc_compare.py +++ b/app/rtc_compare.py @@ -13,11 +13,11 @@ RTC_S1_PRODUCTS_ERROR_ABS_TOLERANCE = 1e-05 LIST_EXCLUDE_COMPARISON = \ - ['//science/CSAR/RTC/metadata/processingInformation/algorithms/ISCEVersion', - '//science/CSAR/RTC/metadata/processingInformation/inputs/auxcalFiles', - '//science/CSAR/RTC/metadata/processingInformation/inputs/configFiles', - '//science/CSAR/RTC/metadata/processingInformation/inputs/demFiles', - '//science/CSAR/RTC/metadata/processingInformation/inputs/orbitFiles'] + ['//science/SENTINEL1/RTC/metadata/processingInformation/algorithms/ISCEVersion', + '//science/SENTINEL1/RTC/metadata/processingInformation/inputs/auxcalFiles', + '//science/SENTINEL1/RTC/metadata/processingInformation/inputs/configFiles', + '//science/SENTINEL1/RTC/metadata/processingInformation/inputs/demFiles', + '//science/SENTINEL1/RTC/metadata/processingInformation/inputs/orbitFiles'] def _get_parser(): @@ -259,9 +259,9 @@ def compare_hdf5_elements(hdf5_obj_1, hdf5_obj_2, str_key, is_attr=False, # convert object reference to the path to which it is pointing # Example: # attribute `REFERENCE_LIST` in - # /science/CSAR/RTC/grids/frequencyA/xCoordinates' + # /science/SENTINEL1/RTC/grids/frequencyA/xCoordinates' # attribute `DIMENSION_LIST` in - # /science/CSAR/RTC/grids/frequencyA/VH + # /science/SENTINEL1/RTC/grids/frequencyA/VH if (len(val_1.shape) >= 1) and ('shape' in dir(val_1[0])): if (isinstance(val_1[0], np.void) or ((len(val_1[0].shape) == 1) and (isinstance(val_1[0][0], h5py.h5r.Reference)))): diff --git a/app/rtc_s1.py b/app/rtc_s1.py index e6cf9e31..3209f4bf 100755 --- a/app/rtc_s1.py +++ b/app/rtc_s1.py @@ -4,6 +4,7 @@ RTC Workflow ''' +import datetime import os import time @@ -18,10 +19,9 @@ from rtc.geogrid import snap_coord from rtc.runconfig import RunConfig -from rtc.mosaic_geobursts import weighted_mosaic -from rtc.core import create_logger -from rtc.h5_prep import save_hdf5_file, create_hdf5_file, \ - save_hdf5_dataset, BASE_DS +from rtc.mosaic_geobursts import compute_weighted_mosaic_raster, compute_weighted_mosaic_raster_single_band +from rtc.core import create_logger, save_as_cog +from rtc.h5_prep import save_hdf5_file, create_hdf5_file, BASE_HDF5_DATASET logger = logging.getLogger('rtc_s1') @@ -67,10 +67,86 @@ def _update_mosaic_boundaries(mosaic_geogrid_dict, geogrid): assert(mosaic_geogrid_dict['epsg'] == geogrid.epsg) +def _separate_pol_channels(multi_band_file, output_file_list, logger, + output_raster_format): + """Save a multi-band raster file as individual single-band files + + Parameters + ---------- + multi_band_file : str + Multi-band raster file + output_file_list : list + Output file list + logger : loggin.Logger + Logger + """ + gdal_ds = gdal.Open(multi_band_file, gdal.GA_ReadOnly) + description = gdal_ds.GetDescription() + projection = gdal_ds.GetProjectionRef() + geotransform = gdal_ds.GetGeoTransform() + metadata = gdal_ds.GetMetadata() + + num_bands = gdal_ds.RasterCount + if num_bands != len(output_file_list): + error_str = (f'ERROR number of output files ({len(output_file_list)})' + f' does not match with the number' + f' of input bursts` bands ({num_bands})') + raise ValueError(error_str) + + for b, output_file in enumerate(output_file_list): + gdal_band = gdal_ds.GetRasterBand(b + 1) + gdal_dtype = gdal_band.DataType + band_image = gdal_band.ReadAsArray() + + # Save the corrected image + driver_out = gdal.GetDriverByName(output_raster_format) + raster_out = driver_out.Create( + output_file, band_image.shape[1], + band_image.shape[0], 1, gdal_dtype) + + raster_out.SetDescription(description) + raster_out.SetProjection(projection) + raster_out.SetGeoTransform(geotransform) + raster_out.SetMetadata(metadata) + + band_out = raster_out.GetRasterBand(1) + band_out.WriteArray(band_image) + band_out.FlushCache() + del band_out + logger.info(f'file saved: {output_file}') + + def _create_raster_obj(output_dir, ds_name, ds_hdf5, dtype, shape, - radar_grid_file_dict, output_obj_list, flag_save_vector_1, + radar_grid_file_dict, output_obj_list, flag_create_raster_obj, extension): - if flag_save_vector_1 is not True: + """Create an ISCE3 raster object (GTiff) for a radar geometry layer. + + Parameters + ---------- + output_dir: str + Output directory + ds_name: str + Dataset (geometry layer) name + ds_hdf5: str + HDF5 dataset name + dtype:: gdal.DataType + GDAL data type + shape: list + Shape of the output raster + radar_grid_file_dict: dict + Dictionary that will hold the name of the output file + referenced by the contents of `ds_hdf5` (dict key) + output_obj_list: list + Mutable list of output raster objects + flag_create_raster_obj: bool + Flag indicating if raster object should be created + + Returns + ------- + raster_obj : isce3.io.Raster + ISCE3 raster object + """ + if flag_create_raster_obj is not True: return None output_file = os.path.join(output_dir, ds_name) + '.' + extension @@ -87,12 +163,12 @@ def _create_raster_obj(output_dir, ds_name, ds_hdf5, dtype, shape, def _add_output_to_output_metadata_dict(flag, key, output_dir, - output_metadata_dict, product_id, extension): + output_metadata_dict, product_prefix, extension): if not flag: return output_image_list = [] output_metadata_dict[key] = \ - [os.path.join(output_dir, f'{product_id}_{key}.{extension}'), + [os.path.join(output_dir, f'{product_prefix}_{key}.{extension}'), output_image_list] @@ -112,7 +188,7 @@ def apply_slc_corrections(burst_in: Sentinel1BurstSlc, # Apply the correction if flag_thermal_correction: - logger.info(f'applying thermal noise correction to burst SLCs') + logger.info(f' applying thermal noise correction to burst SLC') corrected_image = np.abs(arr_slc_from) ** 2 - burst_in.thermal_noise_lut min_backscatter = 0 max_backscatter = None @@ -122,7 +198,7 @@ def apply_slc_corrections(burst_in: Sentinel1BurstSlc, corrected_image=np.abs(arr_slc_from) ** 2 if flag_apply_abs_rad_correction: - logger.info(f'applying absolute radiometric correction to burst SLCs') + logger.info(f' applying absolute radiometric correction to burst SLC') if flag_output_complex: factor_mag = np.sqrt(corrected_image) / np.abs(arr_slc_from) factor_mag[np.isnan(factor_mag)] = 0.0 @@ -147,21 +223,136 @@ def apply_slc_corrections(burst_in: Sentinel1BurstSlc, del band_out -def run(cfg): +def calculate_layover_shadow_mask(burst_in: Sentinel1BurstSlc, + geogrid_in: isce3.product.GeoGridParameters, + path_dem: str, + filename_out: str, + output_raster_format: str, + threshold_rdr2geo: float=1.0e-7, + numiter_rdr2geo: int=25, + extraiter_rdr2geo: int=10, + lines_per_block_rdr2geo: int=1000, + threshold_geo2rdr: float=1.0e-8, + numiter_geo2rdr: int=25, + nlooks_az: int=1, nlooks_rg: int=1): + ''' + Generate the layover shadow mask and geocode the mask + + Parameters: + ----------- + burst_in: Sentinel1BurstSlc + Input burst + geogrid_in: isce3.product.GeoGridParameters + Geogrid to geocode the layover shadow mask in radar grid + path_dem: str + Path to the DEM + filename_out: str + Path to the geocoded layover shadow mask + output_raster_format: str + File format of the layover shadow mask + threshold_rdr2geo: float + Iteration threshold for rdr2geo + numiter_rdr2geo: int + Number of max. iteration for rdr2geo object + extraiter_rdr2geo: int + Extra number of iteration for rdr2geo object + lines_per_block_rdr2geo: int + Lines per block for rdr2geo + threshold_geo2rdr: float + Iteration threshold for geo2rdr + numiter_geo2rdr: int + Number of max. iteration for geo2rdr object + nlooks_az: int + Number of looks in azimuth direction. For the calculation in coarse grid + nlooks_rg: int + Number of looks in range direction. For the calculation in coarse grid + + ''' + + # determine the output filename + str_datetime = burst_in.sensing_start.strftime('%Y%m%d_%H%M%S.%f') + + path_layover_shadow_mask = (f'layover_shadow_mask_{burst_in.burst_id}_' + f'{burst_in.polarization}_{str_datetime}') + + # Run topo to get layover shadow mask + dem_raster = isce3.io.Raster(path_dem) + ellipsoid = isce3.core.Ellipsoid() + + Rdr2Geo = isce3.geometry.Rdr2Geo + + rdr_grid = burst_in.as_isce3_radargrid() + + # when requested, apply mulitilooking on radar grid for the computation in coarse resolution + if nlooks_az > 1 or nlooks_rg > 1: + rdr_grid = rdr_grid.multilook(nlooks_az, nlooks_rg) + + isce3_orbit = burst_in.orbit + grid_doppler = isce3.core.LUT2d() + + rdr2geo_obj = Rdr2Geo(rdr_grid, + isce3_orbit, + ellipsoid, + grid_doppler, + threshold=threshold_rdr2geo, + numiter=numiter_rdr2geo, + extraiter=extraiter_rdr2geo, + lines_per_block=lines_per_block_rdr2geo) + + mask_raster = isce3.io.Raster(path_layover_shadow_mask, rdr_grid.width, + rdr_grid.length, 1, gdal.GDT_Byte, 'MEM') + + rdr2geo_obj.topo(dem_raster, None, None, None, + layover_shadow_raster=mask_raster) + + # geocode the layover shadow mask + geo = isce3.geocode.GeocodeFloat32() + geo.orbit = isce3_orbit + geo.ellipsoid = ellipsoid + geo.doppler = grid_doppler + geo.threshold_geo2rdr = threshold_geo2rdr + geo.numiter_geo2rdr = numiter_geo2rdr + geo.lines_per_block = lines_per_block_rdr2geo + geo.data_interpolator = 'NEAREST' + geo.geogrid(float(geogrid_in.start_x), + float(geogrid_in.start_y), + float(geogrid_in.spacing_x), + float(geogrid_in.spacing_y), + int(geogrid_in.width), + int(geogrid_in.length), + int(geogrid_in.epsg)) + + geocoded_raster = isce3.io.Raster(filename_out, + geogrid_in.width, geogrid_in.length, 1, + gdal.GDT_Byte, output_raster_format) + + geo.geocode(radar_grid=rdr_grid, + input_raster=mask_raster, + output_raster=geocoded_raster, + dem_raster=dem_raster, + output_mode=isce3.geocode.GeocodeOutputMode.INTERP) + + +def run(cfg: RunConfig): ''' Run geocode burst workflow with user-defined args stored in dictionary runconfig `cfg` Parameters --------- - cfg: dict - Dictionary with user runconfig options + cfg: RunConfig + RunConfig object with user runconfig options ''' # Start tracking processing time t_start = time.time() time_stamp = str(float(time.time())) - logger.info("Starting geocode burst") + logger.info("Starting the RTC-S1 Science Application Software (SAS)") + + # primary executable + processing_type = cfg.groups.product_group.processing_type + product_version_float = cfg.groups.product_group.product_version + product_version = f'{product_version_float:.1f}' # unpack processing parameters processing_namespace = cfg.groups.processing @@ -174,33 +365,73 @@ def run(cfg): processing_namespace.apply_absolute_radiometric_correction # read product path group / output format - product_id = cfg.groups.product_path_group.product_id + product_id = cfg.groups.product_group.product_id if product_id is None: product_id = 'rtc_product' + product_prefix = f'{product_id}_v{product_version}' scratch_path = os.path.join( - cfg.groups.product_path_group.scratch_path, f'temp_{time_stamp}') - output_dir = cfg.groups.product_path_group.output_dir - flag_mosaic = cfg.groups.product_path_group.mosaic_bursts + cfg.groups.product_group.scratch_path, f'temp_{time_stamp}') + output_dir = cfg.groups.product_group.output_dir - output_format = cfg.groups.product_path_group.output_format + # RTC-S1 imagery + save_bursts = cfg.groups.product_group.save_bursts + save_mosaics = cfg.groups.product_group.save_mosaics - flag_hdf5 = (output_format == 'HDF5' or output_format == 'NETCDF') + if not save_bursts and not save_mosaics: + err_msg = (f"ERROR either `save_bursts` or `save_mosaics` needs to be" + " set") + raise ValueError(err_msg) - if output_format == 'NETCDF': + output_imagery_format = \ + cfg.groups.product_group.output_imagery_format + output_imagery_compression = \ + cfg.groups.product_group.output_imagery_compression + output_imagery_nbits = \ + cfg.groups.product_group.output_imagery_nbits + + logger.info(f'Identification:') + logger.info(f' product ID: {product_id}') + logger.info(f' processing type: {processing_type}') + logger.info(f' product version: {product_version}') + logger.info(f' product prefix: {product_prefix}') + logger.info(f'Processing parameters:') + logger.info(f' apply RTC: {flag_apply_rtc}') + logger.info(f' apply thermal noise correction:' + f' {flag_apply_thermal_noise_correction}') + logger.info(f' apply absolute radiometric correction:' + f' {flag_apply_abs_rad_correction}') + logger.info(f' scratch dir: {scratch_path}') + logger.info(f' output dir: {output_dir}') + logger.info(f' save bursts: {save_bursts}') + logger.info(f' save mosaics: {save_mosaics}') + logger.info(f' output imagery format: {output_imagery_format}') + logger.info(f' output imagery compression: {output_imagery_compression}') + logger.info(f' output imagery nbits: {output_imagery_nbits}') + + save_imagery_as_hdf5 = (output_imagery_format == 'HDF5' or + output_imagery_format == 'NETCDF') + save_secondary_layers_as_hdf5 = \ + cfg.groups.product_group.save_secondary_layers_as_hdf5 + + save_metadata = (cfg.groups.product_group.save_metadata or + save_imagery_as_hdf5 or + save_secondary_layers_as_hdf5) + + if output_imagery_format == 'NETCDF': hdf5_file_extension = 'nc' else: hdf5_file_extension = 'h5' - if flag_hdf5: + if save_imagery_as_hdf5 or output_imagery_format == 'COG': output_raster_format = 'GTiff' else: - output_raster_format = output_format + output_raster_format = output_imagery_format if output_raster_format == 'GTiff': - extension = 'tif' + imagery_extension = 'tif' else: - extension = 'bin' + imagery_extension = 'bin' # unpack geocode run parameters geocode_namespace = cfg.groups.processing.geocoding @@ -215,22 +446,34 @@ def run(cfg): abs_cal_factor = geocode_namespace.abs_rad_cal clip_max = geocode_namespace.clip_max clip_min = geocode_namespace.clip_min - # geogrids = geocode_namespace.geogrids flag_upsample_radar_grid = geocode_namespace.upsample_radargrid - flag_save_incidence_angle = geocode_namespace.save_incidence_angle - flag_save_local_inc_angle = geocode_namespace.save_local_inc_angle - flag_save_projection_angle = geocode_namespace.save_projection_angle - flag_save_rtc_anf_psi = geocode_namespace.save_rtc_anf_psi - flag_save_range_slope = \ + save_incidence_angle = geocode_namespace.save_incidence_angle + save_local_inc_angle = geocode_namespace.save_local_inc_angle + save_projection_angle = geocode_namespace.save_projection_angle + save_rtc_anf_psi = geocode_namespace.save_rtc_anf_psi + save_range_slope = \ geocode_namespace.save_range_slope - flag_save_nlooks = geocode_namespace.save_nlooks - flag_save_rtc_anf = geocode_namespace.save_rtc_anf - flag_save_dem = geocode_namespace.save_dem + save_nlooks = geocode_namespace.save_nlooks + + + + + + # TODO remove the lines below: + if save_mosaics: + save_nlooks = True + + + + + save_rtc_anf = geocode_namespace.save_rtc_anf + save_dem = geocode_namespace.save_dem + save_layover_shadow_mask = geocode_namespace.save_layover_shadow_mask - flag_call_radar_grid = (flag_save_incidence_angle or - flag_save_local_inc_angle or flag_save_projection_angle or - flag_save_rtc_anf_psi or flag_save_dem or - flag_save_range_slope) + flag_call_radar_grid = (save_incidence_angle or + save_local_inc_angle or save_projection_angle or + save_rtc_anf_psi or save_dem or + save_range_slope) # unpack RTC run parameters rtc_namespace = cfg.groups.processing.rtc @@ -259,32 +502,41 @@ def run(cfg): # Common initializations dem_raster = isce3.io.Raster(cfg.dem) - epsg = dem_raster.get_epsg() - proj = isce3.core.make_projection(epsg) - ellipsoid = proj.ellipsoid + ellipsoid = isce3.core.Ellipsoid() zero_doppler = isce3.core.LUT2d() threshold = cfg.geo2rdr_params.threshold maxiter = cfg.geo2rdr_params.numiter exponent = 1 if (flag_apply_thermal_noise_correction or flag_apply_abs_rad_correction) else 2 - # output mosaics - geo_filename = f'{output_dir}/'f'{product_id}.{extension}' + # output mosaics variables + geo_filename = f'{output_dir}/'f'{product_prefix}.{imagery_extension}' output_imagery_list = [] output_file_list = [] output_metadata_dict = {} - if flag_hdf5: + # output dir (imagery mosaics) + if save_imagery_as_hdf5: output_dir_mosaic_raster = scratch_path else: output_dir_mosaic_raster = output_dir + # output dir (secondary layers mosaics) + if save_secondary_layers_as_hdf5: + output_dir_sec_mosaic_raster = scratch_path + else: + output_dir_sec_mosaic_raster = output_dir + + _add_output_to_output_metadata_dict( + save_layover_shadow_mask, 'layover_shadow_mask', + output_dir_sec_mosaic_raster, + output_metadata_dict, product_prefix, imagery_extension) _add_output_to_output_metadata_dict( - flag_save_nlooks, 'nlooks', output_dir_mosaic_raster, - output_metadata_dict, product_id, extension) + save_nlooks, 'nlooks', output_dir_sec_mosaic_raster, + output_metadata_dict, product_prefix, imagery_extension) _add_output_to_output_metadata_dict( - flag_save_rtc_anf, 'rtc', output_dir_mosaic_raster, - output_metadata_dict, product_id, extension) + save_rtc_anf, 'rtc', output_dir_sec_mosaic_raster, + output_metadata_dict, product_prefix, imagery_extension) mosaic_geogrid_dict = {} temp_files_list = [] @@ -294,11 +546,11 @@ def run(cfg): vrt_options_mosaic = gdal.BuildVRTOptions(separate=True) n_bursts = len(cfg.bursts.items()) - print('number of bursts to process:', n_bursts) + print('Number of bursts to process:', n_bursts) hdf5_obj = None output_hdf5_file = os.path.join(output_dir, - f'{product_id}.{hdf5_file_extension}') + f'{product_prefix}.{hdf5_file_extension}') # iterate over sub-burts for burst_index, (burst_id, burst_pol_dict) in enumerate(cfg.bursts.items()): @@ -306,21 +558,25 @@ def run(cfg): # start burst processing t_burst_start = time.time() - logger.info(f'processing burst: {burst_id} ({burst_index+1}/' + logger.info(f'Processing burst: {burst_id} ({burst_index+1}/' f'{n_bursts})') - pols = list(burst_pol_dict.keys()) - burst = burst_pol_dict[pols[0]] + pol_list = list(burst_pol_dict.keys()) + burst = burst_pol_dict[pol_list[0]] - flag_bursts_files_are_temporary = \ - flag_hdf5 or (flag_mosaic and not n_bursts == 1) + flag_bursts_files_are_temporary = (not save_bursts or + save_imagery_as_hdf5) + flag_bursts_secondary_files_are_temporary = ( + not save_bursts or save_secondary_layers_as_hdf5) burst_scratch_path = f'{scratch_path}/{burst_id}/' os.makedirs(burst_scratch_path, exist_ok=True) - if flag_bursts_files_are_temporary: + if not save_bursts: + # burst files are saved in scratch dir bursts_output_dir = burst_scratch_path else: + # burst files (individual or HDF5) are saved in burst_id dir bursts_output_dir = os.path.join(output_dir, burst_id) os.makedirs(bursts_output_dir, exist_ok=True) @@ -335,7 +591,7 @@ def run(cfg): # update mosaic boundaries _update_mosaic_boundaries(mosaic_geogrid_dict, geogrid) - logger.info(f'reading burst SLCs') + logger.info(f' reading burst SLCs') radar_grid = burst.as_isce3_radargrid() # native_doppler = burst.doppler.lut2d orbit = burst.orbit @@ -347,12 +603,12 @@ def run(cfg): mosaic_geogrid_dict['lookside'] = radar_grid.lookside input_file_list = [] - pol_list = list(burst_pol_dict.keys()) + for pol, burst_pol in burst_pol_dict.items(): temp_slc_path = \ f'{burst_scratch_path}/rslc_{pol}.vrt' temp_slc_corrected_path = ( - f'{burst_scratch_path}/rslc_{pol}_corrected.{extension}') + f'{burst_scratch_path}/rslc_{pol}_corrected.{imagery_extension}') burst_pol.slc_to_vrt_file(temp_slc_path) if (flag_apply_thermal_noise_correction or @@ -366,10 +622,11 @@ def run(cfg): flag_apply_thermal_noise_correction, flag_apply_abs_rad_correction=True) input_burst_filename = temp_slc_corrected_path + temp_files_list.append(temp_slc_corrected_path) else: input_burst_filename = temp_slc_path - temp_files_list.append(input_burst_filename) + temp_files_list.append(temp_slc_path) input_file_list.append(input_burst_filename) # create multi-band VRT @@ -382,18 +639,12 @@ def run(cfg): rdr_burst_raster = isce3.io.Raster(temp_vrt_path) temp_files_list.append(temp_vrt_path) - # Generate output geocoded burst raster - if flag_bursts_files_are_temporary: - # files are temporary - geo_burst_filename = \ - f'{burst_scratch_path}/{product_id}.{extension}' - temp_files_list.append(geo_burst_filename) - else: - os.makedirs(f'{output_dir}/{burst_id}', exist_ok=True) - geo_burst_filename = \ - f'{output_dir}/{burst_id}/{product_id}.{extension}' - output_file_list.append(geo_burst_filename) - + # At this point, burst imagery files are always temporary + geo_burst_filename = \ + f'{burst_scratch_path}/{product_prefix}.{imagery_extension}' + temp_files_list.append(geo_burst_filename) + + # Generate output geocoded burst raster geo_burst_raster = isce3.io.Raster( geo_burst_filename, geogrid.width, geogrid.length, @@ -428,10 +679,10 @@ def run(cfg): geogrid.spacing_x, geogrid.spacing_y, geogrid.width, geogrid.length, geogrid.epsg) - if flag_save_nlooks: - nlooks_file = (f'{bursts_output_dir}/{product_id}' - f'_nlooks.{extension}') - if flag_bursts_files_are_temporary: + if save_nlooks: + nlooks_file = (f'{bursts_output_dir}/{product_prefix}' + f'_nlooks.{imagery_extension}') + if flag_bursts_secondary_files_are_temporary: temp_files_list.append(nlooks_file) else: output_file_list.append(nlooks_file) @@ -442,10 +693,10 @@ def run(cfg): nlooks_file = None out_geo_nlooks_obj = None - if flag_save_rtc_anf: - rtc_anf_file = (f'{bursts_output_dir}/{product_id}' - f'_rtc_anf.{extension}') - if flag_bursts_files_are_temporary: + if save_rtc_anf: + rtc_anf_file = (f'{bursts_output_dir}/{product_prefix}' + f'_rtc_anf.{imagery_extension}') + if flag_bursts_secondary_files_are_temporary: temp_files_list.append(rtc_anf_file) else: output_file_list.append(rtc_anf_file) @@ -475,80 +726,165 @@ def run(cfg): sub_swaths.set_valid_samples_array(1, valid_samples_sub_swath) # geocode - geo_obj.geocode(radar_grid=radar_grid, - input_raster=rdr_burst_raster, - output_raster=geo_burst_raster, - dem_raster=dem_raster, - output_mode=geocode_algorithm, - geogrid_upsampling=geogrid_upsampling, - flag_apply_rtc=flag_apply_rtc, - input_terrain_radiometry=input_terrain_radiometry, - output_terrain_radiometry=output_terrain_radiometry, - exponent=exponent, - rtc_min_value_db=rtc_min_value_db, - rtc_upsampling=rtc_upsampling, - rtc_algorithm=rtc_algorithm, - abs_cal_factor=abs_cal_factor, - flag_upsample_radar_grid=flag_upsample_radar_grid, - clip_min = clip_min, - clip_max = clip_max, - # radargrid_nlooks=radar_grid_nlooks, - # out_off_diag_terms=out_off_diag_terms_obj, - out_geo_nlooks=out_geo_nlooks_obj, - out_geo_rtc=out_geo_rtc_obj, - # out_geo_dem=out_geo_dem_obj, - input_rtc=None, - output_rtc=None, - dem_interp_method=dem_interp_method_enum, - memory_mode=memory_mode, - sub_swaths=sub_swaths) + flag_error_sub_swaths = False + try: + geo_obj.geocode(radar_grid=radar_grid, + input_raster=rdr_burst_raster, + output_raster=geo_burst_raster, + dem_raster=dem_raster, + output_mode=geocode_algorithm, + geogrid_upsampling=geogrid_upsampling, + flag_apply_rtc=flag_apply_rtc, + input_terrain_radiometry=input_terrain_radiometry, + output_terrain_radiometry=output_terrain_radiometry, + exponent=exponent, + rtc_min_value_db=rtc_min_value_db, + rtc_upsampling=rtc_upsampling, + rtc_algorithm=rtc_algorithm, + abs_cal_factor=abs_cal_factor, + flag_upsample_radar_grid=flag_upsample_radar_grid, + clip_min = clip_min, + clip_max = clip_max, + # out_off_diag_terms=out_off_diag_terms_obj, + out_geo_nlooks=out_geo_nlooks_obj, + out_geo_rtc=out_geo_rtc_obj, + input_rtc=None, + output_rtc=None, + dem_interp_method=dem_interp_method_enum, + memory_mode=memory_mode, + sub_swaths=sub_swaths) + except TypeError: + flag_error_sub_swaths = True + logger.warning('WARNING there was an error executing geocode().' + ' Retrying it with less parameters') + + # geocode (without sub_swaths) + geo_obj.geocode(radar_grid=radar_grid, + input_raster=rdr_burst_raster, + output_raster=geo_burst_raster, + dem_raster=dem_raster, + output_mode=geocode_algorithm, + geogrid_upsampling=geogrid_upsampling, + flag_apply_rtc=flag_apply_rtc, + input_terrain_radiometry=input_terrain_radiometry, + output_terrain_radiometry=output_terrain_radiometry, + exponent=exponent, + rtc_min_value_db=rtc_min_value_db, + rtc_upsampling=rtc_upsampling, + rtc_algorithm=rtc_algorithm, + abs_cal_factor=abs_cal_factor, + flag_upsample_radar_grid=flag_upsample_radar_grid, + clip_min = clip_min, + clip_max = clip_max, + # out_off_diag_terms=out_off_diag_terms_obj, + out_geo_nlooks=out_geo_nlooks_obj, + out_geo_rtc=out_geo_rtc_obj, + input_rtc=None, + output_rtc=None, + dem_interp_method=dem_interp_method_enum, + memory_mode=memory_mode) + + if flag_error_sub_swaths: + logger.warning('WARNING the sub-swath masking is not available' + ' from this ISCE3 version. The sub-swath masking' + ' was disabled.') + + # Calculate layover shadow mask when requested + if save_layover_shadow_mask: + layover_shadow_mask_file = (f'{bursts_output_dir}/{product_prefix}' + f'_layover_shadow_mask.{imagery_extension}') + calculate_layover_shadow_mask(burst, + geogrid, + cfg.dem, + layover_shadow_mask_file, + output_raster_format, + threshold_rdr2geo=cfg.rdr2geo_params.threshold, + numiter_rdr2geo=cfg.rdr2geo_params.numiter, + threshold_geo2rdr=cfg.geo2rdr_params.threshold, + numiter_geo2rdr=cfg.geo2rdr_params.numiter) + + if flag_bursts_secondary_files_are_temporary: + temp_files_list.append(layover_shadow_mask_file) + else: + output_file_list.append(layover_shadow_mask_file) + logger.info(f'file saved: {layover_shadow_mask_file}') + output_metadata_dict['layover_shadow_mask'][1].append( + layover_shadow_mask_file) + + else: + layover_shadow_mask_file = None del geo_burst_raster - if not flag_bursts_files_are_temporary: - logger.info(f'file saved: {geo_burst_filename}') + + # Output imagery list contains multi-band files that + # will be used for mosaicking output_imagery_list.append(geo_burst_filename) - if flag_save_nlooks: + # If burst imagery is not temporary, separate polarization channels + if not flag_bursts_files_are_temporary: + output_burst_imagery_list = [] + for pol in pol_list: + geo_burst_pol_filename = \ + os.path.join(output_dir, burst_id, + f'{product_prefix}_{pol}.' + + f'{imagery_extension}') + output_burst_imagery_list.append(geo_burst_pol_filename) + + _separate_pol_channels(geo_burst_filename, + output_burst_imagery_list, + logger, output_raster_format) + + output_file_list += output_burst_imagery_list + + if save_nlooks: del out_geo_nlooks_obj - if not flag_bursts_files_are_temporary: + if not flag_bursts_secondary_files_are_temporary: logger.info(f'file saved: {nlooks_file}') output_metadata_dict['nlooks'][1].append(nlooks_file) - if flag_save_rtc_anf: + if save_rtc_anf: del out_geo_rtc_obj - if not flag_bursts_files_are_temporary: + if not flag_bursts_secondary_files_are_temporary: logger.info(f'file saved: {rtc_anf_file}') output_metadata_dict['rtc'][1].append(rtc_anf_file) radar_grid_file_dict = {} - if flag_call_radar_grid and not flag_mosaic: + if flag_call_radar_grid and save_bursts: get_radar_grid( - geogrid, dem_interp_method_enum, product_id, - bursts_output_dir, extension, flag_save_incidence_angle, - flag_save_local_inc_angle, flag_save_projection_angle, - flag_save_rtc_anf_psi, - flag_save_range_slope, flag_save_dem, + geogrid, dem_interp_method_enum, product_prefix, + bursts_output_dir, imagery_extension, save_incidence_angle, + save_local_inc_angle, save_projection_angle, + save_rtc_anf_psi, + save_range_slope, save_dem, dem_raster, radar_grid_file_dict, mosaic_geogrid_dict, orbit, - verbose = not flag_bursts_files_are_temporary) - if flag_hdf5: + verbose = not flag_bursts_secondary_files_are_temporary) + if flag_bursts_secondary_files_are_temporary: # files are temporary temp_files_list += list(radar_grid_file_dict.values()) else: output_file_list += list(radar_grid_file_dict.values()) - if flag_hdf5 and not flag_mosaic: + # Create burst HDF5 + if ((save_imagery_as_hdf5 or save_metadata) and save_bursts): hdf5_file_output_dir = os.path.join(output_dir, burst_id) os.makedirs(hdf5_file_output_dir, exist_ok=True) - output_hdf5_file = os.path.join( - hdf5_file_output_dir, f'{product_id}.{hdf5_file_extension}') - hdf5_obj = create_hdf5_file(output_hdf5_file, orbit, burst, cfg) + output_hdf5_file_burst = os.path.join( + hdf5_file_output_dir, f'{product_prefix}.{hdf5_file_extension}') + hdf5_obj = create_hdf5_file(output_hdf5_file_burst, orbit, burst, cfg) save_hdf5_file( - hdf5_obj, output_hdf5_file, flag_apply_rtc, - clip_max, clip_min, output_radiometry_str, output_file_list, + hdf5_obj, output_hdf5_file_burst, flag_apply_rtc, + clip_max, clip_min, output_radiometry_str, geogrid, pol_list, geo_burst_filename, nlooks_file, - rtc_anf_file, radar_grid_file_dict) - elif flag_hdf5 and flag_mosaic and burst_index == 0: + rtc_anf_file, layover_shadow_mask_file, + radar_grid_file_dict, + save_imagery = save_imagery_as_hdf5, + save_secondary_layers = save_secondary_layers_as_hdf5) + output_file_list.append(output_hdf5_file_burst) + + # Create mosaic HDF5 + if ((save_imagery_as_hdf5 or save_metadata) and save_mosaics + and burst_index == 0): hdf5_obj = create_hdf5_file(output_hdf5_file, orbit, burst, cfg) t_burst_end = time.time() @@ -558,61 +894,85 @@ def run(cfg): # end burst processing # =========================================================== - if flag_call_radar_grid and flag_mosaic: + if flag_call_radar_grid and save_mosaics: radar_grid_file_dict = {} - if flag_hdf5: + if save_secondary_layers_as_hdf5: radar_grid_output_dir = scratch_path else: radar_grid_output_dir = output_dir - get_radar_grid(cfg.geogrid, dem_interp_method_enum, product_id, - radar_grid_output_dir, extension, flag_save_incidence_angle, - flag_save_local_inc_angle, flag_save_projection_angle, - flag_save_rtc_anf_psi, - flag_save_range_slope, flag_save_dem, + get_radar_grid(cfg.geogrid, dem_interp_method_enum, product_prefix, + radar_grid_output_dir, imagery_extension, save_incidence_angle, + save_local_inc_angle, save_projection_angle, + save_rtc_anf_psi, + save_range_slope, save_dem, dem_raster, radar_grid_file_dict, mosaic_geogrid_dict, - orbit, verbose = not flag_hdf5) - if flag_hdf5: + orbit, verbose = not save_imagery_as_hdf5) + if save_imagery_as_hdf5: # files are temporary temp_files_list += list(radar_grid_file_dict.values()) else: output_file_list += list(radar_grid_file_dict.values()) - if flag_mosaic: - # mosaic sub-bursts - geo_filename = f'{output_dir_mosaic_raster}/{product_id}.{extension}' - logger.info(f'mosaicking file: {geo_filename}') + if save_mosaics: + + # Mosaic sub-bursts imagery + logger.info(f'mosaicking files:') + output_imagery_filename_list = [] + for pol in pol_list: + geo_pol_filename = \ + (f'{output_dir_mosaic_raster}/{product_prefix}_{pol}.' + f'{imagery_extension}') + logger.info(f' {geo_pol_filename}') + output_imagery_filename_list.append(geo_pol_filename) nlooks_list = output_metadata_dict['nlooks'][1] - weighted_mosaic(output_imagery_list, nlooks_list, - geo_filename, cfg.geogrid, verbose=False) + compute_weighted_mosaic_raster_single_band( + output_imagery_list, nlooks_list, + output_imagery_filename_list, cfg.geogrid, verbose=False) - if flag_hdf5: - temp_files_list.append(geo_filename) + if save_imagery_as_hdf5: + temp_files_list += output_imagery_filename_list else: - output_file_list.append(geo_filename) + output_file_list += output_imagery_filename_list - # mosaic other bands + # Mosaic other bands for key in output_metadata_dict.keys(): output_file, input_files = output_metadata_dict[key] logger.info(f'mosaicking file: {output_file}') - weighted_mosaic(input_files, nlooks_list, output_file, + compute_weighted_mosaic_raster(input_files, nlooks_list, output_file, cfg.geogrid, verbose=False) - if flag_hdf5: + + + + + # TODO: Remove nlooks exception below + if (save_secondary_layers_as_hdf5 or + (key == 'nlooks' and not save_nlooks)): temp_files_list.append(output_file) else: output_file_list.append(output_file) - if flag_hdf5: - if flag_save_nlooks: + + + + + # Save HDF5 + if save_imagery_as_hdf5 or save_metadata: + if save_nlooks: nlooks_mosaic_file = output_metadata_dict['nlooks'][0] else: nlooks_mosaic_file = None - if flag_save_rtc_anf: + if save_rtc_anf: rtc_anf_mosaic_file = output_metadata_dict['rtc'][0] else: rtc_anf_mosaic_file = None + if save_layover_shadow_mask: + layover_shadow_mask_file = output_metadata_dict[ + 'layover_shadow_mask'][0] + else: + layover_shadow_mask_file = None # Update metadata datasets that depend on all bursts sensing_start = None @@ -620,21 +980,16 @@ def run(cfg): for burst_id, burst_pol_dict in cfg.bursts.items(): pols = list(burst_pol_dict.keys()) burst = burst_pol_dict[pols[0]] - print('this burst:') - if sensing_start is not None: - print(' ', sensing_start.strftime('%Y-%m-%dT%H:%M:%S.%f')) - if sensing_stop is not None: - print(' ', sensing_stop.strftime('%Y-%m-%dT%H:%M:%S.%f')) + if (sensing_start is None or burst.sensing_start < sensing_start): sensing_start = burst.sensing_start - print('updated sensing start') + if sensing_stop is None or burst.sensing_stop > sensing_stop: sensing_stop = burst.sensing_stop - print('updated sensing stop') - sensing_start_ds = f'{BASE_DS}/identification/zeroDopplerStartTime' - sensing_end_ds = f'{BASE_DS}/identification/zeroDopplerEndTime' + sensing_start_ds = f'{BASE_HDF5_DATASET}/identification/zeroDopplerStartTime' + sensing_end_ds = f'{BASE_HDF5_DATASET}/identification/zeroDopplerEndTime' del hdf5_obj[sensing_start_ds] del hdf5_obj[sensing_end_ds] hdf5_obj[sensing_start_ds] = \ @@ -642,11 +997,24 @@ def run(cfg): hdf5_obj[sensing_end_ds] = \ sensing_stop.strftime('%Y-%m-%dT%H:%M:%S.%f') - save_hdf5_file(hdf5_obj, output_hdf5_file, flag_apply_rtc, - clip_max, clip_min, output_radiometry_str, - output_file_list, cfg.geogrid, pol_list, - geo_filename, nlooks_mosaic_file, - rtc_anf_mosaic_file, radar_grid_file_dict) + save_hdf5_file( + hdf5_obj, output_hdf5_file, flag_apply_rtc, + clip_max, clip_min, output_radiometry_str, + cfg.geogrid, pol_list, geo_filename, nlooks_mosaic_file, + rtc_anf_mosaic_file, layover_shadow_mask_file, + radar_grid_file_dict, save_imagery = save_imagery_as_hdf5, + save_secondary_layers = save_secondary_layers_as_hdf5) + output_file_list.append(output_hdf5_file) + + if output_imagery_format == 'COG': + logger.info(f'Saving files as Cloud-Optimized GeoTIFFs (COGs)') + for filename in output_file_list: + if not filename.endswith('.tif'): + continue + logger.info(f' processing file: {filename}') + save_as_cog(filename, scratch_path, logger, + compression=output_imagery_compression, + nbits=output_imagery_nbits) logger.info('removing temporary files:') for filename in temp_files_list: @@ -664,11 +1032,11 @@ def run(cfg): -def get_radar_grid(geogrid, dem_interp_method_enum, product_id, - output_dir, extension, flag_save_incidence_angle, - flag_save_local_inc_angle, flag_save_projection_angle, - flag_save_rtc_anf_psi, - flag_save_range_slope, flag_save_dem, dem_raster, +def get_radar_grid(geogrid, dem_interp_method_enum, product_prefix, + output_dir, extension, save_incidence_angle, + save_local_inc_angle, save_projection_angle, + save_rtc_anf_psi, + save_range_slope, save_dem, dem_raster, radar_grid_file_dict, mosaic_geogrid_dict, orbit, verbose = True): output_obj_list = [] @@ -676,31 +1044,31 @@ def get_radar_grid(geogrid, dem_interp_method_enum, product_id, shape = [layers_nbands, geogrid.length, geogrid.width] incidence_angle_raster = _create_raster_obj( - output_dir, f'{product_id}_incidence_angle', + output_dir, f'{product_prefix}_incidence_angle', 'incidenceAngle', gdal.GDT_Float32, shape, radar_grid_file_dict, - output_obj_list, flag_save_incidence_angle, extension) + output_obj_list, save_incidence_angle, extension) local_incidence_angle_raster = _create_raster_obj( - output_dir, f'{product_id}_local_incidence_angle', + output_dir, f'{product_prefix}_local_incidence_angle', 'localIncidenceAngle', gdal.GDT_Float32, shape, - radar_grid_file_dict, output_obj_list, flag_save_local_inc_angle, + radar_grid_file_dict, output_obj_list, save_local_inc_angle, extension) projection_angle_raster = _create_raster_obj( - output_dir, f'{product_id}_projection_angle', + output_dir, f'{product_prefix}_projection_angle', 'projectionAngle', gdal.GDT_Float32, shape, radar_grid_file_dict, - output_obj_list, flag_save_projection_angle, extension) + output_obj_list, save_projection_angle, extension) rtc_anf_psi_raster = _create_raster_obj( - output_dir, f'{product_id}_rtc_anf_psi', + output_dir, f'{product_prefix}_rtc_anf_psi', 'areaNormalizationFactorPsi', gdal.GDT_Float32, shape, radar_grid_file_dict, output_obj_list, - flag_save_rtc_anf_psi, extension) + save_rtc_anf_psi, extension) range_slope_raster = _create_raster_obj( - output_dir, f'{product_id}_range_slope', + output_dir, f'{product_prefix}_range_slope', 'rangeSlope', gdal.GDT_Float32, shape, radar_grid_file_dict, - output_obj_list, flag_save_range_slope, extension) + output_obj_list, save_range_slope, extension) interpolated_dem_raster = _create_raster_obj( - output_dir, f'{product_id}_interpolated_dem', + output_dir, f'{product_prefix}_interpolated_dem', 'interpolatedDem', gdal.GDT_Float32, shape, radar_grid_file_dict, - output_obj_list, flag_save_dem, extension) + output_obj_list, save_dem, extension) # TODO review this (Doppler)!!! # native_doppler = burst.doppler.lut2d @@ -709,6 +1077,13 @@ def get_radar_grid(geogrid, dem_interp_method_enum, product_id, grid_doppler = isce3.core.LUT2d() grid_doppler.bounds_error = False + # TODO: update code below + # Computation of range slope is not merged to ISCE yet + kwargs_get_radar_grid = {} + if range_slope_raster: + kwargs_get_radar_grid['directional_slope_angle_raster'] = \ + range_slope_raster + # call get_radar_grid() isce3.geogrid.get_radar_grid(mosaic_geogrid_dict['lookside'], mosaic_geogrid_dict['wavelength'], @@ -725,11 +1100,10 @@ def get_radar_grid(geogrid, dem_interp_method_enum, product_id, projection_angle_raster, simulated_radar_brightness_raster = rtc_anf_psi_raster, - directional_slope_angle_raster = - range_slope_raster, interpolated_dem_raster = interpolated_dem_raster, - dem_interp_method=dem_interp_method_enum) + dem_interp_method=dem_interp_method_enum, + **kwargs_get_radar_grid) # Flush data for obj in output_obj_list: diff --git a/build_docker_image.sh b/build_docker_image.sh index 4e7e94d1..f20b0ccb 100755 --- a/build_docker_image.sh +++ b/build_docker_image.sh @@ -2,7 +2,7 @@ REPO=opera IMAGE=rtc -TAG=interface_0.1 +TAG=bet_0.2 echo "IMAGE is $REPO/$IMAGE:$TAG" @@ -19,7 +19,7 @@ docker build --rm --force-rm --network host -t $REPO/$IMAGE:$TAG -f Docker/Docke #docker run --rm -u "$(id -u):$(id -g)" -v "$PWD:/mnt" -w /mnt -it --network host "${IMAGE}:$t" pytest /mnt/tests/ # create image tar -docker save $REPO/$IMAGE:$TAG > Docker/dockerimg_rtc_interface_0.1.tar +docker save $REPO/$IMAGE:$TAG > Docker/dockerimg_rtc_beta_0.2.tar # remove image docker image rm $REPO/$IMAGE:$TAG diff --git a/setup.py b/setup.py index 96088adb..136745db 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,8 @@ classifiers = ['Programming Language :: Python',], scripts = ['app/rtc_s1.py'], install_requires = ['argparse', 'numpy', 'yamale', - 'scipy', 'pytest', 'requests'], + 'scipy', 'pytest', 'requests', + 'pyproj', 'shapely'], url = 'https://github.com/opera-adt/RTC', author = ('Gustavo H. X. Shiroma', 'Seongsu Jeong'), author_email = ('gustavo.h.shiroma@jpl.nasa.gov, seongsu.jeong@jpl.nasa.gov'), diff --git a/src/rtc/core.py b/src/rtc/core.py index 03607dbe..84059f9c 100644 --- a/src/rtc/core.py +++ b/src/rtc/core.py @@ -1,7 +1,12 @@ #!/usr/bin/env python +import os import sys +import shutil import logging +import tempfile + +from osgeo import gdal logger = logging.getLogger('rtc_s1') @@ -47,6 +52,104 @@ def flush(self): self.buffer = '' + +def save_as_cog(filename, scratch_dir = '.', logger = None, + flag_compress=True, ovr_resamp_algorithm=None, + compression='DEFLATE', nbits=None): + """Save (overwrite) a GeoTIFF file as a cloud-optimized GeoTIFF. + + Parameters + ---------- + filename: str + GeoTIFF to be saved as a cloud-optimized GeoTIFF + scratch_dir: str (optional) + Temporary Directory + ovr_resamp_algorithm: str (optional) + Resampling algorithm for overviews. + Options: "AVERAGE", "AVERAGE_MAGPHASE", "RMS", "BILINEAR", + "CUBIC", "CUBICSPLINE", "GAUSS", "LANCZOS", "MODE", + "NEAREST", or "NONE". Defaults to "NEAREST", if integer, and + "CUBICSPLINE", otherwise. + compression: str (optional) + Compression type. + Optional: "NONE", "LZW", "JPEG", "DEFLATE", "ZSTD", "WEBP", + "LERC", "LERC_DEFLATE", "LERC_ZSTD", "LZMA" + + """ + if logger is None: + logger = logging.getLogger('proteus') + + logger.info(' COG step 1: add overviews') + gdal_ds = gdal.Open(filename, gdal.GA_Update) + gdal_dtype = gdal_ds.GetRasterBand(1).DataType + dtype_name = gdal.GetDataTypeName(gdal_dtype).lower() + + overviews_list = [4, 16, 64, 128] + + is_integer = 'byte' in dtype_name or 'int' in dtype_name + if ovr_resamp_algorithm is None and is_integer: + ovr_resamp_algorithm = 'NEAREST' + elif ovr_resamp_algorithm is None: + ovr_resamp_algorithm = 'CUBICSPLINE' + + gdal_ds.BuildOverviews(ovr_resamp_algorithm, overviews_list, + gdal.TermProgress_nocb) + + del gdal_ds # close the dataset (Python object and pointers) + external_overview_file = filename + '.ovr' + if os.path.isfile(external_overview_file): + os.remove(external_overview_file) + + logger.info(' COG step 2: save as COG') + temp_file = tempfile.NamedTemporaryFile( + dir=scratch_dir, suffix='.tif').name + + # Blocks of 512 x 512 => 256 KiB (UInt8) or 1MiB (Float32) + tile_size = 512 + gdal_translate_options = ['BIGTIFF=IF_SAFER', + 'MAX_Z_ERROR=0', + 'TILED=YES', + f'BLOCKXSIZE={tile_size}', + f'BLOCKYSIZE={tile_size}', + 'COPY_SRC_OVERVIEWS=YES'] + + if compression: + gdal_translate_options += [f'COMPRESS={compression}'] + + if is_integer: + gdal_translate_options += ['PREDICTOR=2'] + else: + gdal_translate_options += ['PREDICTOR=3'] + + if nbits is not None: + gdal_translate_options += [f'NBITS={nbits}'] + + # suppress type casting errors + gdal.SetConfigOption('CPL_LOG', '/dev/null') + + gdal.Translate(temp_file, filename, + creationOptions=gdal_translate_options) + + shutil.move(temp_file, filename) + + logger.info(' COG step 3: validate') + try: + from rtc.extern.validate_cloud_optimized_geotiff import main as validate_cog + except ModuleNotFoundError: + logger.info('WARNING could not import module validate_cloud_optimized_geotiff') + return + + argv = ['--full-check=yes', filename] + validate_cog_ret = validate_cog(argv) + if validate_cog_ret == 0: + logger.info(f' file "{filename}" is a valid cloud optimized' + ' GeoTIFF') + else: + logger.warning(f' file "{filename}" is NOT a valid cloud' + f' optimized GeoTIFF!') + + + def create_logger(log_file, full_log_formatting=None): """Create logger object for a log file diff --git a/src/rtc/defaults/rtc_s1.yaml b/src/rtc/defaults/rtc_s1.yaml index 65e21874..0a5148f1 100644 --- a/src/rtc/defaults/rtc_s1.yaml +++ b/src/rtc/defaults/rtc_s1.yaml @@ -1,147 +1,194 @@ runconfig: - name: rtc_s1_workflow_default - - groups: - pge_name_group: - pge_name: RTC_S1_PGE - - input_file_group: - # Required. List of SAFE files (min=1) - safe_file_path: - # Required. List of orbit (EOF) files (min=1) - orbit_file_path: - # Required. The unique burst ID to process - burst_id: - - dynamic_ancillary_file_group: - # Digital elevation model - dem_file: - - product_path_group: - # Directory where PGE will place results - product_path: - # Directory where SAS writes temporary data - scratch_path: - - # If option `mosaic_bursts` is not set, output files are saved to: - # {output_dir}/{burst_id}/{product_id}{suffix}.{ext} - # If option `mosaic_bursts` is set, output files are saved to: - # {output_dir}/{product_id}{suffix}.{ext} - # If the field `product_id`` is left empty, the prefix "rtc_product" - # will be used instead. - # `suffix` is only used when there are multiple output files. - # `ext` is determined by geocoding_options.output_format. - output_dir: - product_id: - mosaic_bursts: False - - # Format of output file - output_format: NETCDF - - - primary_executable: - product_type: RTC_S1 - - processing: - polarization: - geo2rdr: - threshold: 1.0e-8 - numiter: 25 - - dem_interpolation_method: biquintic - - # Apply absolute radiometric correction - apply_absolute_radiometric_correction: True + name: rtc_s1_workflow_default + + groups: + pge_name_group: + pge_name: RTC_S1_PGE + + input_file_group: + # Required. List of SAFE files (min=1) + safe_file_path: + # Required. List of orbit (EOF) files (min=1) + orbit_file_path: + # Optional. Burst ID to process (empty for all bursts) + burst_id: + + dynamic_ancillary_file_group: + # Digital elevation model + dem_file: + + # Digital elevation model source description + dem_description: + + static_ancillary_file_group: + + # burst database sqlite file + burst_database_file: + + product_group: + + processing_type: UNDEFINED + + product_version: 0.2 + + # Directory where PGE will place results + product_path: + # Directory where SAS writes temporary data + scratch_path: + + # If option `save_bursts` is set, output bursts are saved to: + # {output_dir}/{burst_id}/{product_id}_v{product_version}{suffix}.{ext} + # If option `save_mosaics` is set, output mosaics are saved to: + # {output_dir}/{product_id}_v{product_version}{suffix}.{ext} + # If the field `product_id`` is left empty, the prefix "rtc_product" + # will be used instead. + # `suffix` is only used when there are multiple output files. + # `ext` is determined by geocoding_options.output_imagery_format. + output_dir: + product_id: + + # RTC-S1 imagery + save_bursts: True + save_mosaics: False + output_imagery_format: COG + output_imagery_compression: ZSTD + output_imagery_nbits: 16 + + # Optional. Save secondary layers (e.g., inc. angle) within + # the HDF5 file + save_secondary_layers_as_hdf5: True + + # Save RTC-S1 metadata in the HDF5 format + # Optional for `output_imagery_format` equal to 'ENVI', 'GTiff', or + # 'COG', and enabled by default for `output_imagery_format` equal + # to 'HDF5' or 'NETCDF' or `save_secondary_layers_as_hdf5` is True + save_metadata: True + + primary_executable: + product_type: RTC_S1 + + processing: + + # Check if ancillary inputs cover entirely the output product + check_ancillary_inputs_coverage: True + + # Polarization channels to process. + polarization: + + # Options to run geo2rdr + geo2rdr: + threshold: 1.0e-7 + numiter: 50 + + # Options to run rdr2geo + rdr2geo: + threshold: 1.0e-7 + numiter: 25 + + # DEM interpolation method + dem_interpolation_method: biquintic + + # Apply absolute radiometric correction + apply_absolute_radiometric_correction: True - # Apply thermal noise correction - apply_thermal_noise_correction: True + # Apply thermal noise correction + apply_thermal_noise_correction: True + + # OPTIONAL - Apply RTC + apply_rtc: True + + # Apply bistatic delay correction + apply_bistatic_delay_correction: True - # OPTIONAL - Apply RTC - apply_rtc: True + # Apply dry tropospheric delay correction + apply_dry_tropospheric_delay_correction: True - # OPTIONAL - to control behavior of RTC module - # (only applicable if geocode.apply_rtc is True) - rtc: - # OPTIONAL - Choices: - # "gamma0" (default) - # "sigma0" - output_type: gamma0 - - # OPTIONAL - Choices: - # "bilinear_distribution" (default) - # "area_projection" - algorithm_type: area_projection - - # OPTIONAL - Choices: - # "beta0" (default) - # "sigma0" - input_terrain_radiometry: beta0 - - # OPTIONAL - Minimum RTC area factor in dB - rtc_min_value_db: - - # RTC DEM upsampling - dem_upsampling: 1 - - # Geocoding options - geocoding: - - # OPTIONAL - Algorithm type, area projection or - # interpolation: sinc, bilinear, bicubic, nearest, and biquintic - algorithm_type: area_projection + # OPTIONAL - to control behavior of RTC module + # (only applicable if geocode.apply_rtc is True) + rtc: + # OPTIONAL - Choices: + # "gamma0" (default) + # "sigma0" + output_type: gamma0 + + # OPTIONAL - Choices: + # "bilinear_distribution" (default) + # "area_projection" + algorithm_type: area_projection + + # OPTIONAL - Choices: + # "beta0" (default) + # "sigma0" + input_terrain_radiometry: beta0 + + # OPTIONAL - Minimum RTC area factor in dB + rtc_min_value_db: -40 + + # RTC DEM upsampling + dem_upsampling: 2 + + # Geocoding options + geocoding: + + # OPTIONAL - Algorithm type, area projection or + # interpolation: sinc, bilinear, bicubic, nearest, and biquintic + algorithm_type: area_projection - # OPTIONAL - Choices: "single_block", "geogrid", "geogrid_and_radargrid", and "auto" (default) - memory_mode: + # OPTIONAL - Choices: "single_block", "geogrid", "geogrid_and_radargrid", and "auto" (default) + memory_mode: - # OPTIONAL - Processing upsampling factor applied to input geogrid - geogrid_upsampling: 1 + # OPTIONAL - Processing upsampling factor applied to input geogrid + geogrid_upsampling: 1 - # Save the incidence angle - save_incidence_angle: False + # Save the incidence angle + save_incidence_angle: False - # Save the local-incidence angle - save_local_inc_angle: False + # Save the local-incidence angle + save_local_inc_angle: False - # Save the projection angle - save_projection_angle: False + # Save the projection angle + save_projection_angle: False - # Save the RTC area normalization factor (ANF) computed with - # the projection angle method - save_rtc_anf_psi: False + # Save the RTC area normalization factor (ANF) computed with + # the projection angle method + save_rtc_anf_psi: False - # Save the range slope angle - save_range_slope: False + # Save the range slope angle + save_range_slope: False - # Save the number of looks used to generate the RTC product - save_nlooks: True + # Save the number of looks used to generate the RTC product + save_nlooks: False - # Save the RTC area normalization factor (ANF) used to generate - # the RTC product - save_rtc_anf: True - - # Save the interpolated DEM used to generate the RTC product - save_dem: False - - # OPTIONAL - Absolute radiometric correction - abs_rad_cal: 1 - - # OPTIONAL - Clip values above threshold - clip_max: - - # OPTIONAL - Clip values below threshold - clip_min: - - # Double SLC sampling in the range direction - upsample_radargrid: False - - output_epsg: - x_posting: 30 - y_posting: 30 - x_snap: 30 - y_snap: 30 - top_left: - x: - y: - bottom_right: - x: - y: + # Save the RTC area normalization factor (ANF) used to generate + # the RTC product + save_rtc_anf: False + + # Save the interpolated DEM used to generate the RTC product + save_dem: False + + # Save layover shadow mask + save_layover_shadow_mask: True + + # OPTIONAL - Absolute radiometric correction + abs_rad_cal: 1 + + # OPTIONAL - Clip values above threshold + clip_max: + + # OPTIONAL - Clip values below threshold + clip_min: + + # Double SLC sampling in the range direction + upsample_radargrid: False + + output_epsg: + x_posting: 30 + y_posting: 30 + x_snap: 30 + y_snap: 30 + top_left: + x: + y: + bottom_right: + x: + y: diff --git a/src/rtc/extern/__init__.py b/src/rtc/extern/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/rtc/extern/validate_cloud_optimized_geotiff.py b/src/rtc/extern/validate_cloud_optimized_geotiff.py new file mode 100644 index 00000000..c8855f06 --- /dev/null +++ b/src/rtc/extern/validate_cloud_optimized_geotiff.py @@ -0,0 +1,407 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# ***************************************************************************** +# $Id$ +# +# Project: GDAL +# Purpose: Validate Cloud Optimized GeoTIFF file structure +# Author: Even Rouault, +# +# ***************************************************************************** +# Copyright (c) 2017, Even Rouault +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the "Software"), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# ***************************************************************************** + +import os.path +import struct +import sys +from osgeo import gdal + + +def Usage(): + print('Usage: validate_cloud_optimized_geotiff.py [-q] [--full-check=yes/no/auto] test.tif') + print('') + print('Options:') + print('-q: quiet mode') + print('--full-check=yes/no/auto: check tile/strip leader/trailer bytes. auto=yes for local files, and no for remote files') + return 1 + + +class ValidateCloudOptimizedGeoTIFFException(Exception): + pass + + +def full_check_band(f, band_name, band, errors, + block_order_row_major, + block_leader_size_as_uint4, + block_trailer_last_4_bytes_repeated, + mask_interleaved_with_imagery): + + block_size = band.GetBlockSize() + mask_band = None + if mask_interleaved_with_imagery: + mask_band = band.GetMaskBand() + mask_block_size = mask_band.GetBlockSize() + if block_size != mask_block_size: + errors += [band_name + ': mask block size is different from its imagery band'] + mask_band = None + + yblocks = (band.YSize + block_size[1] - 1) // block_size[1] + xblocks = (band.XSize + block_size[0] - 1) // block_size[0] + last_offset = 0 + for y in range(yblocks): + for x in range(xblocks): + + offset = band.GetMetadataItem('BLOCK_OFFSET_%d_%d' % (x, y), 'TIFF') + offset = int(offset) if offset is not None else 0 + bytecount = band.GetMetadataItem('BLOCK_SIZE_%d_%d' % (x, y), 'TIFF') + bytecount = int(bytecount) if bytecount is not None else 0 + + if offset > 0: + if block_order_row_major and offset < last_offset: + errors += [band_name + + ': offset of block (%d, %d) is smaller than previous block' % (x, y)] + + if block_leader_size_as_uint4: + gdal.VSIFSeekL(f, offset - 4, 0) + leader_size = struct.unpack('= 4: + gdal.VSIFSeekL(f, offset + bytecount - 4, 0) + last_bytes = gdal.VSIFReadL(8, 1, f) + if last_bytes[0:4] != last_bytes[4:8]: + errors += [band_name + + ': for block (%d, %d), trailer bytes are invalid' % (x, y)] + + if mask_band: + offset_mask = mask_band.GetMetadataItem('BLOCK_OFFSET_%d_%d' % (x, y), 'TIFF') + offset_mask = int(offset_mask) if offset_mask is not None else 0 + if offset > 0 and offset_mask > 0: + #bytecount_mask = int(mask_band.GetMetadataItem('BLOCK_SIZE_%d_%d' % (x,y), 'TIFF')) + expected_offset_mask = offset + bytecount + \ + (4 if block_leader_size_as_uint4 else 0) + \ + (4 if block_trailer_last_4_bytes_repeated else 0) + if offset_mask != expected_offset_mask: + errors += ['Mask of ' + band_name + ': for block (%d, %d), offset is %d, whereas %d was expected' % ( + x, y, offset_mask, expected_offset_mask)] + elif offset == 0 and offset_mask > 0: + if block_order_row_major and offset_mask < last_offset: + errors += ['Mask of ' + band_name + + ': offset of block (%d, %d) is smaller than previous block' % (x, y)] + + offset = offset_mask + + last_offset = offset + + +def validate(ds, check_tiled=True, full_check=False): + """Check if a file is a (Geo)TIFF with cloud optimized compatible structure. + + Args: + ds: GDAL Dataset for the file to inspect. + check_tiled: Set to False to ignore missing tiling. + full_check: Set to TRUe to check tile/strip leader/trailer bytes. Might be slow on remote files + + Returns: + A tuple, whose first element is an array of error messages + (empty if there is no error), and the second element, a dictionary + with the structure of the GeoTIFF file. + + Raises: + ValidateCloudOptimizedGeoTIFFException: Unable to open the file or the + file is not a Tiff. + """ + + if int(gdal.VersionInfo('VERSION_NUM')) < 2020000: + raise ValidateCloudOptimizedGeoTIFFException( + 'GDAL 2.2 or above required') + + unicode_type = type(''.encode('utf-8').decode('utf-8')) + if isinstance(ds, (str, unicode_type)): + gdal.PushErrorHandler() + ds = gdal.Open(ds) + gdal.PopErrorHandler() + if ds is None: + raise ValidateCloudOptimizedGeoTIFFException( + 'Invalid file : %s' % gdal.GetLastErrorMsg()) + if ds.GetDriver().ShortName != 'GTiff': + raise ValidateCloudOptimizedGeoTIFFException( + 'The file is not a GeoTIFF') + + details = {} + errors = [] + warnings = [] + filename = ds.GetDescription() + main_band = ds.GetRasterBand(1) + ovr_count = main_band.GetOverviewCount() + filelist = ds.GetFileList() + if filelist is not None and filename + '.ovr' in filelist: + errors += [ + 'Overviews found in external .ovr file. They should be internal'] + + if main_band.XSize > 512 or main_band.YSize > 512: + if check_tiled: + block_size = main_band.GetBlockSize() + if block_size[0] == main_band.XSize and block_size[0] > 1024: + errors += [ + 'The file is greater than 512xH or Wx512, but is not tiled'] + + if ovr_count == 0: + warnings += [ + 'The file is greater than 512xH or Wx512, it is recommended ' + 'to include internal overviews'] + + ifd_offset = int(main_band.GetMetadataItem('IFD_OFFSET', 'TIFF')) + ifd_offsets = [ifd_offset] + + block_order_row_major = False + block_leader_size_as_uint4 = False + block_trailer_last_4_bytes_repeated = False + mask_interleaved_with_imagery = False + + if ifd_offset not in (8, 16): + + # Check if there is GDAL hidden structural metadata + f = gdal.VSIFOpenL(filename, 'rb') + if not f: + raise ValidateCloudOptimizedGeoTIFFException("Cannot open file") + signature = struct.unpack('B' * 4, gdal.VSIFReadL(4, 1, f)) + bigtiff = signature in ((0x49, 0x49, 0x2B, 0x00), (0x4D, 0x4D, 0x00, 0x2B)) + if bigtiff: + expected_ifd_pos = 16 + else: + expected_ifd_pos = 8 + gdal.VSIFSeekL(f, expected_ifd_pos, 0) + pattern = "GDAL_STRUCTURAL_METADATA_SIZE=%06d bytes\n" % 0 + got = gdal.VSIFReadL(len(pattern), 1, f).decode('LATIN1') + if len(got) == len(pattern) and got.startswith('GDAL_STRUCTURAL_METADATA_SIZE='): + size = int(got[len('GDAL_STRUCTURAL_METADATA_SIZE='):][0:6]) + extra_md = gdal.VSIFReadL(size, 1, f).decode('LATIN1') + block_order_row_major = 'BLOCK_ORDER=ROW_MAJOR' in extra_md + block_leader_size_as_uint4 = 'BLOCK_LEADER=SIZE_AS_UINT4' in extra_md + block_trailer_last_4_bytes_repeated = 'BLOCK_TRAILER=LAST_4_BYTES_REPEATED' in extra_md + mask_interleaved_with_imagery = 'MASK_INTERLEAVED_WITH_IMAGERY=YES' in extra_md + if 'KNOWN_INCOMPATIBLE_EDITION=YES' in extra_md: + errors += ["KNOWN_INCOMPATIBLE_EDITION=YES is declared in the file"] + expected_ifd_pos += len(pattern) + size + expected_ifd_pos += expected_ifd_pos % 2 # IFD offset starts on a 2-byte boundary + gdal.VSIFCloseL(f) + + if expected_ifd_pos != ifd_offsets[0]: + errors += [ + 'The offset of the main IFD should be %d. It is %d instead' % (expected_ifd_pos, ifd_offsets[0])] + + details['ifd_offsets'] = {} + details['ifd_offsets']['main'] = ifd_offset + + for i in range(ovr_count): + # Check that overviews are by descending sizes + ovr_band = ds.GetRasterBand(1).GetOverview(i) + if i == 0: + if (ovr_band.XSize > main_band.XSize or + ovr_band.YSize > main_band.YSize): + errors += [ + 'First overview has larger dimension than main band'] + else: + prev_ovr_band = ds.GetRasterBand(1).GetOverview(i - 1) + if (ovr_band.XSize > prev_ovr_band.XSize or + ovr_band.YSize > prev_ovr_band.YSize): + errors += [ + 'Overview of index %d has larger dimension than ' + 'overview of index %d' % (i, i - 1)] + + if check_tiled: + block_size = ovr_band.GetBlockSize() + if block_size[0] == ovr_band.XSize and block_size[0] > 1024: + errors += [ + 'Overview of index %d is not tiled' % i] + + # Check that the IFD of descending overviews are sorted by increasing + # offsets + ifd_offset = int(ovr_band.GetMetadataItem('IFD_OFFSET', 'TIFF')) + ifd_offsets.append(ifd_offset) + details['ifd_offsets']['overview_%d' % i] = ifd_offset + if ifd_offsets[-1] < ifd_offsets[-2]: + if i == 0: + errors += [ + 'The offset of the IFD for overview of index %d is %d, ' + 'whereas it should be greater than the one of the main ' + 'image, which is at byte %d' % + (i, ifd_offsets[-1], ifd_offsets[-2])] + else: + errors += [ + 'The offset of the IFD for overview of index %d is %d, ' + 'whereas it should be greater than the one of index %d, ' + 'which is at byte %d' % + (i, ifd_offsets[-1], i - 1, ifd_offsets[-2])] + + # Check that the imagery starts by the smallest overview and ends with + # the main resolution dataset + + def get_block_offset(band): + blockxsize, blockysize = band.GetBlockSize() + for y in range(int((band.YSize + blockysize - 1) / blockysize)): + for x in range(int((band.XSize + blockxsize - 1) / blockxsize)): + block_offset = band.GetMetadataItem('BLOCK_OFFSET_%d_%d' % (x, y), 'TIFF') + if block_offset: + return int(block_offset) + return 0 + + block_offset = get_block_offset(main_band) + data_offsets = [block_offset] + details['data_offsets'] = {} + details['data_offsets']['main'] = block_offset + for i in range(ovr_count): + ovr_band = ds.GetRasterBand(1).GetOverview(i) + block_offset = get_block_offset(ovr_band) + data_offsets.append(block_offset) + details['data_offsets']['overview_%d' % i] = block_offset + + if data_offsets[-1] != 0 and data_offsets[-1] < ifd_offsets[-1]: + if ovr_count > 0: + errors += [ + 'The offset of the first block of the smallest overview ' + 'should be after its IFD'] + else: + errors += [ + 'The offset of the first block of the image should ' + 'be after its IFD'] + for i in range(len(data_offsets) - 2, 0, -1): + if data_offsets[i] != 0 and data_offsets[i] < data_offsets[i + 1]: + errors += [ + 'The offset of the first block of overview of index %d should ' + 'be after the one of the overview of index %d' % + (i - 1, i)] + if len(data_offsets) >= 2 and data_offsets[0] != 0 and data_offsets[0] < data_offsets[1]: + errors += [ + 'The offset of the first block of the main resolution image ' + 'should be after the one of the overview of index %d' % + (ovr_count - 1)] + + if full_check and (block_order_row_major or block_leader_size_as_uint4 or + block_trailer_last_4_bytes_repeated or + mask_interleaved_with_imagery): + f = gdal.VSIFOpenL(filename, 'rb') + if not f: + raise ValidateCloudOptimizedGeoTIFFException("Cannot open file") + + full_check_band(f, 'Main resolution image', main_band, errors, + block_order_row_major, + block_leader_size_as_uint4, + block_trailer_last_4_bytes_repeated, + mask_interleaved_with_imagery) + if main_band.GetMaskFlags() == gdal.GMF_PER_DATASET and \ + (filename + '.msk') not in ds.GetFileList(): + full_check_band(f, 'Mask band of main resolution image', + main_band.GetMaskBand(), errors, + block_order_row_major, + block_leader_size_as_uint4, + block_trailer_last_4_bytes_repeated, False) + for i in range(ovr_count): + ovr_band = ds.GetRasterBand(1).GetOverview(i) + full_check_band(f, 'Overview %d' % i, ovr_band, errors, + block_order_row_major, + block_leader_size_as_uint4, + block_trailer_last_4_bytes_repeated, + mask_interleaved_with_imagery) + if ovr_band.GetMaskFlags() == gdal.GMF_PER_DATASET and \ + (filename + '.msk') not in ds.GetFileList(): + full_check_band(f, 'Mask band of overview %d' % i, + ovr_band.GetMaskBand(), errors, + block_order_row_major, + block_leader_size_as_uint4, + block_trailer_last_4_bytes_repeated, False) + gdal.VSIFCloseL(f) + + return warnings, errors, details + + +def main(argv): + """Return 0 in case of success, 1 for failure.""" + + i = 1 + filename = None + quiet = False + full_check = None + while i < len(argv): + if argv[i] == '-q': + quiet = True + elif argv[i] == '--full-check=yes': + full_check = True + elif argv[i] == '--full-check=no': + full_check = False + elif argv[i] == '--full-check=auto': + full_check = None + elif argv[i][0] == '-': + return Usage() + elif filename is None: + filename = argv[i] + else: + return Usage() + + i += 1 + + if filename is None: + return Usage() + + if full_check is None: + full_check = filename.startswith('/vsimem/') or os.path.exists(filename) + + try: + ret = 0 + warnings, errors, details = validate(filename, full_check=full_check) + if warnings: + if not quiet: + print('The following warnings were found:') + for warning in warnings: + print(' - ' + warning) + print('') + if errors: + if not quiet: + print('%s is NOT a valid cloud optimized GeoTIFF.' % filename) + print('The following errors were found:') + for error in errors: + print(' - ' + error) + print('') + ret = 1 + else: + if not quiet: + print('%s is a valid cloud optimized GeoTIFF' % filename) + + if not quiet and not warnings and not errors: + headers_size = min(details['data_offsets'][k] for k in details['data_offsets']) + if headers_size == 0: + headers_size = gdal.VSIStatL(filename).size + print('\nThe size of all IFD headers is %d bytes' % headers_size) + except ValidateCloudOptimizedGeoTIFFException as e: + if not quiet: + print('%s is NOT a valid cloud optimized GeoTIFF : %s' % + (filename, str(e))) + ret = 1 + + return ret + + +if __name__ == '__main__': + sys.exit(main(sys.argv)) diff --git a/src/rtc/geogrid.py b/src/rtc/geogrid.py index ac63793c..bade0b31 100644 --- a/src/rtc/geogrid.py +++ b/src/rtc/geogrid.py @@ -8,10 +8,13 @@ from nisar.workflows.geogrid import _grid_size import isce3 +from rtc.helpers import burst_bboxes_from_db + logger = logging.getLogger('rtc_s1') -def assign_check_geogrid(geogrid, x_start=None, y_start=None, - x_end=None, y_end=None): + +def assign_check_geogrid(geogrid, xmin=None, ymax=None, + xmax=None, ymin=None): ''' Initialize geogrid with user defined parameters. Check the validity of user-defined parameters @@ -20,13 +23,13 @@ def assign_check_geogrid(geogrid, x_start=None, y_start=None, ---------- geogrid: isce3.product.GeoGridParameters ISCE3 object defining the geogrid - x_start: float + xmin: float Geogrid top-left X coordinate - y_start: float + ymax: float Geogrid top-left Y coordinate - x_end: float + xmax: float Geogrid bottom-right X coordinate - y_end: float + ymin: float Geogrid bottom-right Y coordinate Returns @@ -36,29 +39,29 @@ def assign_check_geogrid(geogrid, x_start=None, y_start=None, ''' # Check assigned input coordinates and initialize geogrid accordingly - if x_start is not None: + if xmin is not None: current_end_x = geogrid.start_x + geogrid.spacing_x * geogrid.width - geogrid.start_x = x_start - geogrid.width = int(np.ceil((current_end_x - x_start) / + geogrid.start_x = xmin + geogrid.width = int(np.ceil((current_end_x - xmin) / geogrid.spacing_x)) # Restore geogrid end point if provided by the user - if x_end is not None: - geogrid.width = int(np.ceil((x_end - geogrid.start_x) / + if xmax is not None: + geogrid.width = int(np.ceil((xmax - geogrid.start_x) / geogrid.spacing_x)) - if y_start is not None: + if ymax is not None: current_end_y = geogrid.start_y + geogrid.spacing_y * geogrid.length - geogrid.start_y = y_start - geogrid.length = int(np.ceil((current_end_y - y_start) / + geogrid.start_y = ymax + geogrid.length = int(np.ceil((current_end_y - ymax) / geogrid.spacing_y)) - if y_end is not None: - geogrid.length = int(np.ceil((y_end - geogrid.start_y) / + if ymin is not None: + geogrid.length = int(np.ceil((ymin - geogrid.start_y) / geogrid.spacing_y)) return geogrid -def intersect_geogrid(geogrid, x_start=None, y_start=None, - x_end=None, y_end=None): +def intersect_geogrid(geogrid, xmin=None, ymax=None, + xmax=None, ymin=None): ''' Return intersected geogrid with user defined parameters. @@ -66,13 +69,13 @@ def intersect_geogrid(geogrid, x_start=None, y_start=None, ---------- geogrid: isce3.product.GeoGridParameters ISCE3 object defining the geogrid - x_start: float + xmin: float Geogrid top-left X coordinate - y_start: float + ymax: float Geogrid top-left Y coordinate - x_end: float + xmax: float Geogrid bottom-right X coordinate - y_end: float + ymin: float Geogrid bottom-right Y coordinate Returns @@ -81,31 +84,31 @@ def intersect_geogrid(geogrid, x_start=None, y_start=None, ISCE3 geogrid ''' - if x_start is not None and x_start > geogrid.start_x: + if xmin is not None and xmin > geogrid.start_x: current_end_x = geogrid.start_x + geogrid.spacing_x * geogrid.width - geogrid.start_x = x_start - geogrid.width = int(np.ceil((current_end_x - x_start) / + geogrid.start_x = xmin + geogrid.width = int(np.ceil((current_end_x - xmin) / geogrid.spacing_x)) - if (x_end is not None and - (x_end < geogrid.start_x + geogrid.width * geogrid.spacing_x)): - geogrid.width = int(np.ceil((x_end - geogrid.start_x) / + if (xmax is not None and + (xmax < geogrid.start_x + geogrid.width * geogrid.spacing_x)): + geogrid.width = int(np.ceil((xmax - geogrid.start_x) / geogrid.spacing_x)) - if y_start is not None and y_start < geogrid.start_y: + if ymax is not None and ymax < geogrid.start_y: current_end_y = geogrid.start_y + geogrid.spacing_y * geogrid.length - geogrid.start_y = y_start - geogrid.length = int(np.ceil((current_end_y - y_start) / + geogrid.start_y = ymax + geogrid.length = int(np.ceil((current_end_y - ymax) / geogrid.spacing_y)) - if (y_end is not None and - (y_end > geogrid.start_y + geogrid.length * geogrid.spacing_y)): - geogrid.length = int(np.ceil((y_end - geogrid.start_y) / + if (ymin is not None and + (ymin > geogrid.start_y + geogrid.length * geogrid.spacing_y)): + geogrid.length = int(np.ceil((ymin - geogrid.start_y) / geogrid.spacing_y)) return geogrid -def check_geogrid_endpoints(geogrid, x_end=None, y_end=None): +def check_geogrid_endpoints(geogrid, xmax=None, ymin=None): ''' Check validity of geogrid end points @@ -113,25 +116,26 @@ def check_geogrid_endpoints(geogrid, x_end=None, y_end=None): ----------- geogrid: isce3.product.GeoGridParameters ISCE3 object defining the geogrid - x_end: float + xmax: float Geogrid bottom right X coordinate - y_end: float + ymin: float Geogrid bottom right Y coordinates Returns ------- - x_end: float + xmax: float Verified geogrid bottom-right X coordinate - y_end: float + ymin: float Verified geogrid bottom-right Y coordinate ''' end_pt = lambda start, sz, spacing: start + spacing * sz - if x_end is None: - x_end = end_pt(geogrid.start_x, geogrid.spacing_x, geogrid.width) - if y_end is None: - y_end = end_pt(geogrid.start_y, geogrid.spacing_y, geogrid.length) - return x_end, y_end + if xmax is None: + xmax = end_pt(geogrid.start_x, geogrid.spacing_x, geogrid.width) + if ymin is None: + ymin = end_pt(geogrid.start_y, geogrid.spacing_y, geogrid.length) + + return xmax, ymin def check_snap_values(x_snap, y_snap, x_spacing, y_spacing): @@ -213,16 +217,16 @@ def snap_geogrid(geogrid, x_snap, y_snap): geogrid: isce3.product.GeoGridParameters ISCE3 object containing the snapped geogrid ''' - x_end = geogrid.start_x + geogrid.width * geogrid.spacing_x - y_end = geogrid.start_y + geogrid.length * geogrid.spacing_y + xmax = geogrid.start_x + geogrid.width * geogrid.spacing_x + ymin = geogrid.start_y + geogrid.length * geogrid.spacing_y if x_snap is not None: geogrid.start_x = snap_coord(geogrid.start_x, x_snap, np.floor) - end_x = snap_coord(x_end, x_snap, np.ceil) + end_x = snap_coord(xmax, x_snap, np.ceil) geogrid.width = _grid_size(end_x, geogrid.start_x, geogrid.spacing_x) if y_snap is not None: geogrid.start_y = snap_coord(geogrid.start_y, y_snap, np.ceil) - end_y = snap_coord(y_end, y_snap, np.floor) + end_y = snap_coord(ymin, y_snap, np.floor) geogrid.length = _grid_size(end_y, geogrid.start_y, geogrid.spacing_y) return geogrid @@ -264,33 +268,219 @@ def get_point_epsg(lat, lon): raise ValueError(err_str) + + + + + + +def generate_geogrids_from_db(bursts, geo_dict, dem, burst_db_file): + ''' + Compute frame and bursts geogrids + + Parameters + ---------- + bursts: list[s1reader.s1_burst_slc.Sentinel1BurstSlc] + List of S-1 burst SLC objects + geo_dict: dict + Dictionary containing runconfig processing.geocoding + parameters + dem_file: str + Dem file + burst_db_file : str + Location of burst database sqlite file + + Returns + ------- + geogrid_all: isce3.product.GeoGridParameters + Frame geogrid + geogrids_dict: dict + Dict containing bursts geogrids indexed by burst_id + ''' + + # TODO use `dem_raster` to update `epsg`, if not provided in the runconfig + # dem_raster = isce3.io.Raster(dem) + + # Unpack values from geocoding dictionary + epsg_runconfig = geo_dict['output_epsg'] + xmin_runconfig = geo_dict['top_left']['x'] + ymax_runconfig = geo_dict['top_left']['y'] + x_spacing = geo_dict['x_posting'] + y_spacing_positive = geo_dict['y_posting'] + xmax_runconfig = geo_dict['bottom_right']['x'] + ymin_runconfig = geo_dict['bottom_right']['y'] + x_snap = geo_dict['x_snap'] + y_snap = geo_dict['y_snap'] + + xmin_all_bursts = np.inf + ymax_all_bursts = -np.inf + xmax_all_bursts = -np.inf + ymin_all_bursts = np.inf + + # Check spacing in X direction + if x_spacing is not None and x_spacing <= 0: + err_str = 'Pixel spacing in X/longitude direction needs to be positive' + err_str += f' (x_spacing: {x_spacing})' + logger.error(err_str) + raise ValueError(err_str) + + # Check spacing in Y direction + if y_spacing_positive is not None and y_spacing_positive <= 0: + err_str = 'Pixel spacing in Y/latitude direction needs to be positive' + err_str += f'(y_spacing: {y_spacing_positive})' + logger.error(err_str) + raise ValueError(err_str) + elif y_spacing_positive is not None: + y_spacing = - y_spacing_positive + else: + y_spacing = None + + geogrids_dict = {} + + # get all burst IDs and their EPSGs + bounding boxes + burst_ids = [b for b in bursts] + epsg_bbox_dict = burst_bboxes_from_db(burst_ids, burst_db_file) + + for burst_id, burst_pol in bursts.items(): + pol_list = list(burst_pol.keys()) + burst = burst_pol[pol_list[0]] + if burst_id in geogrids_dict.keys(): + continue + + radar_grid = burst.as_isce3_radargrid() + orbit = burst.orbit + + # extract EPSG and bbox for current burst from dict + # bottom right = (xmax, ymin) and top left = (xmin, ymax) + epsg, (xmin, ymin, xmax, ymax) = epsg_bbox_dict[burst_id] + + if epsg_runconfig and epsg != epsg_runconfig: + err_str = 'ERROR runconfig and burst-DB EPSG codes do not match:' + err_str += f' runconfig: {epsg_runconfig}.' + err_str += f' burst-DB: {epsg}' + logger.error(err_str) + raise NotImplementedError(err_str) + + if x_spacing is None and epsg == 4326: + x_spacing = -0.00027 + elif x_spacing is None: + x_spacing = -30 + + if y_spacing_positive is None and epsg == 4326: + y_spacing = -0.00027 + elif y_spacing_positive is None: + y_spacing = -30 + + # Initialize geogrid with estimated boundaries + geogrid_burst = isce3.product.bbox_to_geogrid( + radar_grid, orbit, isce3.core.LUT2d(), x_spacing, y_spacing, + epsg) + + width = int(np.ceil((xmax - xmin) / x_spacing)) + length = int(np.ceil((ymin - ymax) / y_spacing)) + + geogrid_burst = isce3.product.GeoGridParameters( + xmin, ymax, x_spacing, y_spacing, width, + length, epsg) + + # intersect burst geogrid with runconfig parameters + geogrid_burst = intersect_geogrid( + geogrid_burst, xmin_runconfig, ymax_runconfig, xmax_runconfig, + ymin_runconfig) + + # Check snap values + check_snap_values(x_snap, y_snap, x_spacing, y_spacing) + + # Snap coordinates + geogrid = snap_geogrid(geogrid_burst, x_snap, y_snap) + + xmin_all_bursts = min([xmin_all_bursts, geogrid.start_x]) + ymax_all_bursts = max([ymax_all_bursts, geogrid.start_y]) + xmax_all_bursts = max([xmax_all_bursts, + geogrid.start_x + geogrid.spacing_x * + geogrid.width]) + ymin_all_bursts = min([ymin_all_bursts, + geogrid.start_y + geogrid.spacing_y * + geogrid.length]) + + geogrids_dict[burst_id] = geogrid + + # make sure that user's runconfig parameters have precedence + if xmin_runconfig is not None: + xmin_all_bursts = xmin_runconfig + if ymax_runconfig is not None: + ymax_all_bursts = ymax_runconfig + if xmax_runconfig is not None: + xmax_all_bursts = xmax_runconfig + if ymin_runconfig is not None: + ymin_all_bursts = ymin_runconfig + + width_all = _grid_size(xmax_all_bursts, xmin_all_bursts, x_spacing) + length_all = _grid_size(ymin_all_bursts, ymax_all_bursts, y_spacing) + geogrid_all = isce3.product.GeoGridParameters( + xmin_all_bursts, ymax_all_bursts, x_spacing, y_spacing, width_all, + length_all, epsg) + + # Check snap values + check_snap_values(x_snap, y_snap, x_spacing, y_spacing) + + # Snap coordinates + geogrid_all = snap_geogrid(geogrid_all, x_snap, y_snap) + return geogrid_all, geogrids_dict + + + + + + def generate_geogrids(bursts, geo_dict, dem): + ''' + Compute frame and bursts geogrids + + Parameters + ---------- + bursts: list[s1reader.s1_burst_slc.Sentinel1BurstSlc] + List of S-1 burst SLC objects + geo_dict: dict + Dictionary containing runconfig processing.geocoding + parameters + dem_file: str + Dem file + + Returns + ------- + geogrid_all: isce3.product.GeoGridParameters + Frame geogrid + geogrids_dict: dict + Dict containing bursts geogrids indexed by burst_id + ''' - dem_raster = isce3.io.Raster(dem) + # TODO use `dem_raster` to update `epsg`, if not provided in the runconfig + # dem_raster = isce3.io.Raster(dem) - # Unpack values from geocoding disctionary + # Unpack values from geocoding dictionary epsg = geo_dict['output_epsg'] - x_start = geo_dict['top_left']['x'] - y_start = geo_dict['top_left']['y'] + xmin = geo_dict['top_left']['x'] + ymax = geo_dict['top_left']['y'] x_spacing = geo_dict['x_posting'] y_spacing_positive = geo_dict['y_posting'] - x_end = geo_dict['bottom_right']['x'] - y_end = geo_dict['bottom_right']['y'] + xmax = geo_dict['bottom_right']['x'] + ymin = geo_dict['bottom_right']['y'] x_snap = geo_dict['x_snap'] y_snap = geo_dict['y_snap'] - x_start_all_bursts = np.inf - y_start_all_bursts = -np.inf - x_end_all_bursts = -np.inf - y_end_all_bursts = np.inf + xmin_all_bursts = np.inf + ymax_all_bursts = -np.inf + xmax_all_bursts = -np.inf + ymin_all_bursts = np.inf # Compute burst EPSG if not assigned in runconfig if epsg is None: y_list = [] x_list = [] for burst_pol in bursts.values(): - pols = list(burst_pol.keys()) - burst = burst_pol[pols[0]] + pol_list = list(burst_pol.keys()) + burst = burst_pol[pol_list[0]] y_list.append(burst.center.y) x_list.append(burst.center.x) y_mean = np.nanmean(y_list) @@ -324,34 +514,34 @@ def generate_geogrids(bursts, geo_dict, dem): geogrids_dict = {} for burst_id, burst_pol in bursts.items(): - pols = list(burst_pol.keys()) - burst = burst_pol[pols[0]] - if burst_id in geogrids_dict: + pol_list = list(burst_pol.keys()) + burst = burst_pol[pol_list[0]] + if burst_id in geogrids_dict.keys(): continue radar_grid = burst.as_isce3_radargrid() orbit = burst.orbit geogrid_burst = None - if len(bursts) > 1 or None in [x_start, y_start, x_end, y_end]: + if len(bursts) > 1 or None in [xmin, ymax, xmax, ymin]: # Initialize geogrid with estimated boundaries geogrid_burst = isce3.product.bbox_to_geogrid( radar_grid, orbit, isce3.core.LUT2d(), x_spacing, y_spacing, epsg) - if len(bursts) == 1 and None in [x_start, y_start, x_end, y_end]: + if len(bursts) == 1 and None in [xmin, ymax, xmax, ymin]: # Check and further initialize geogrid geogrid_burst = assign_check_geogrid( - geogrid_burst, x_start, y_start, x_end, y_end) + geogrid_burst, xmin, ymax, xmax, ymin) geogrid_burst = intersect_geogrid( - geogrid_burst, x_start, y_start, x_end, y_end) + geogrid_burst, xmin, ymax, xmax, ymin) else: # If all the start/end coordinates have been assigned, # initialize the geogrid with them - width = _grid_size(x_end, x_start, x_spacing) - length = _grid_size(y_end, y_start, y_spacing) + width = _grid_size(xmax, xmin, x_spacing) + length = _grid_size(ymin, ymax, y_spacing) geogrid_burst = isce3.product.GeoGridParameters( - x_start, y_start, x_spacing, y_spacing, width, length, + xmin, ymax, x_spacing, y_spacing, width, length, epsg) # Check snap values @@ -359,30 +549,30 @@ def generate_geogrids(bursts, geo_dict, dem): # Snap coordinates geogrid = snap_geogrid(geogrid_burst, x_snap, y_snap) - x_start_all_bursts = min([x_start_all_bursts, geogrid.start_x]) - y_start_all_bursts = max([y_start_all_bursts, geogrid.start_y]) - x_end_all_bursts = max([x_end_all_bursts, + xmin_all_bursts = min([xmin_all_bursts, geogrid.start_x]) + ymax_all_bursts = max([ymax_all_bursts, geogrid.start_y]) + xmax_all_bursts = max([xmax_all_bursts, geogrid.start_x + geogrid.spacing_x * geogrid.width]) - y_end_all_bursts = min([y_end_all_bursts, + ymin_all_bursts = min([ymin_all_bursts, geogrid.start_y + geogrid.spacing_y * geogrid.length]) geogrids_dict[burst_id] = geogrid - if x_start is None: - x_start = x_start_all_bursts - if y_start is None: - y_start = y_start_all_bursts - if x_end is None: - x_end = x_end_all_bursts - if y_end is None: - y_end = y_end_all_bursts - - width = _grid_size(x_end, x_start, x_spacing) - length = _grid_size(y_end, y_start, y_spacing) + if xmin is None: + xmin = xmin_all_bursts + if ymax is None: + ymax = ymax_all_bursts + if xmax is None: + xmax = xmax_all_bursts + if ymin is None: + ymin = ymin_all_bursts + + width = _grid_size(xmax, xmin, x_spacing) + length = _grid_size(ymin, ymax, y_spacing) geogrid_all = isce3.product.GeoGridParameters( - x_start, y_start, x_spacing, y_spacing, width, length, epsg) + xmin, ymax, x_spacing, y_spacing, width, length, epsg) # Check snap values check_snap_values(x_snap, y_snap, x_spacing, y_spacing) diff --git a/src/rtc/h5_prep.py b/src/rtc/h5_prep.py index 72913ebb..b8c15d51 100644 --- a/src/rtc/h5_prep.py +++ b/src/rtc/h5_prep.py @@ -7,23 +7,58 @@ from osgeo import gdal import isce3 +import shapely from s1reader.s1_burst_slc import Sentinel1BurstSlc +from s1reader.version import release_version from rtc.runconfig import RunConfig from nisar.workflows.h5_prep import set_get_geo_info -BASE_DS = f'/science/CSAR' +# Version: BETA v0.2 +__version__ = 0.2 + + +BASE_HDF5_DATASET = f'/science/SENTINEL1' FREQ_GRID_SUB_PATH = 'RTC/grids/frequencyA' -FREQ_GRID_DS = f'{BASE_DS}/{FREQ_GRID_SUB_PATH}' +FREQ_GRID_DS = f'{BASE_HDF5_DATASET}/{FREQ_GRID_SUB_PATH}' logger = logging.getLogger('rtc_s1') +def get_polygon_wkt(burst_in: Sentinel1BurstSlc): + ''' + Get WKT for butst's bounding polygon + It returns "POLYGON" when + there is only one polygon that defines the burst's border + It returns "MULTIPOLYGON" when + there is more than one polygon that defines the burst's border + + Parameters: + ----------- + burst_in: Sentinel1BurstSlc + Input burst + + Return: + _ : str + "POLYGON" or "MULTIPOLYGON" in WKT + as the bounding polygon of the input burst + + ''' + + if len(burst_in.border) ==1: + geometry_polygon = burst_in.border[0] + else: + geometry_polygon = shapely.geometry.MultiPolygon(burst_in.border) + + return geometry_polygon.wkt + def save_hdf5_file(hdf5_obj, output_hdf5_file, flag_apply_rtc, clip_max, - clip_min, output_radiometry_str, output_file_list, + clip_min, output_radiometry_str, geogrid, pol_list, geo_burst_filename, nlooks_file, - rtc_anf_file, radar_grid_file_dict): + rtc_anf_file, layover_shadow_mask_file, + radar_grid_file_dict, + save_imagery=True, save_secondary_layers=True): # save grids metadata h5_ds = os.path.join(FREQ_GRID_DS, 'listOfPolarizations') @@ -39,38 +74,50 @@ def save_hdf5_file(hdf5_obj, output_hdf5_file, flag_apply_rtc, clip_max, del hdf5_obj[h5_ds] dset = hdf5_obj.create_dataset(h5_ds, data=bool(flag_apply_rtc)) + + # save geogrid coordinates yds, xds = set_get_geo_info(hdf5_obj, FREQ_GRID_DS, geogrid) - # save RTC imagery - save_hdf5_dataset(geo_burst_filename, hdf5_obj, FREQ_GRID_DS, - yds, xds, pol_list, - long_name=output_radiometry_str, - units='', - valid_min=clip_min, - valid_max=clip_max) - # save nlooks - if nlooks_file: - save_hdf5_dataset(nlooks_file, hdf5_obj, FREQ_GRID_DS, - yds, xds, 'numberOfLooks', - long_name = 'number of looks', - units = '', - valid_min = 0) - - # save rtc - if rtc_anf_file: - save_hdf5_dataset(rtc_anf_file, hdf5_obj, FREQ_GRID_DS, - yds, xds, 'areaNormalizationFactor', - long_name = 'RTC area factor', - units = '', - valid_min = 0) - - for ds_hdf5, filename in radar_grid_file_dict.items(): - save_hdf5_dataset(filename, hdf5_obj, FREQ_GRID_DS, yds, xds, ds_hdf5, - long_name = '', units = '') + if save_imagery: + # save RTC imagery + save_hdf5_dataset(geo_burst_filename, hdf5_obj, FREQ_GRID_DS, + yds, xds, pol_list, + long_name=output_radiometry_str, + units='', + valid_min=clip_min, + valid_max=clip_max) + + if save_secondary_layers: + # save nlooks + if nlooks_file: + save_hdf5_dataset(nlooks_file, hdf5_obj, FREQ_GRID_DS, + yds, xds, 'numberOfLooks', + long_name = 'number of looks', + units = '', + valid_min = 0) + + # save rtc + if rtc_anf_file: + save_hdf5_dataset(rtc_anf_file, hdf5_obj, FREQ_GRID_DS, + yds, xds, 'areaNormalizationFactor', + long_name = 'RTC area factor', + units = '', + valid_min = 0) + + # save layover shadow mask + if layover_shadow_mask_file: + save_hdf5_dataset(layover_shadow_mask_file, hdf5_obj, FREQ_GRID_DS, + yds, xds, 'layoverShadowMask', + long_name = 'Layover/shadow mask', + units = '', + valid_min = 0) + + for ds_hdf5, filename in radar_grid_file_dict.items(): + save_hdf5_dataset(filename, hdf5_obj, FREQ_GRID_DS, yds, xds, ds_hdf5, + long_name = '', units = '') logger.info(f'file saved: {output_hdf5_file}') - output_file_list.append(output_hdf5_file) def create_hdf5_file(output_hdf5_file, orbit, burst, cfg): @@ -85,7 +132,8 @@ def create_hdf5_file(output_hdf5_file, orbit, burst, cfg): populate_metadata_group(hdf5_obj, burst, cfg) # save orbit - orbit_group = hdf5_obj.require_group(f'{BASE_DS}/RTC/metadata/orbit') + orbit_group = hdf5_obj.require_group( + f'{BASE_HDF5_DATASET}/RTC/metadata/orbit') save_orbit(orbit, orbit_group) return hdf5_obj @@ -102,6 +150,9 @@ def save_orbit(orbit, orbit_group): orbit_group["velocity"].attrs["description"] = np.string_("Velocity vector" " record. This record contains the platform velocity data with" " respect to WGS84 G1762 reference frame") + orbit_group.create_dataset( + 'referenceEpoch', + data=np.string_(orbit.reference_epoch.isoformat())) # Orbit source/type # TODO: Update orbit type: @@ -113,7 +164,7 @@ def save_orbit(orbit, orbit_group): def populate_metadata_group(h5py_obj: h5py.File, burst_in: Sentinel1BurstSlc, cfg_in: RunConfig, - root_path: str = BASE_DS): + root_path: str = BASE_HDF5_DATASET): '''Populate RTC metadata based on Sentinel1BurstSlc and RunConfig Parameters: @@ -127,9 +178,26 @@ def populate_metadata_group(h5py_obj: h5py.File, root_path: str Root path inside the HDF5 object on which the metadata will be placed ''' + + # orbit files orbit_files = [os.path.basename(f) for f in cfg_in.orbit_path] + + # L1 SLC granules l1_slc_granules = [os.path.basename(f) for f in cfg_in.safe_files] - dem_files = [os.path.basename(cfg_in.dem)] + + # processing type + processing_type = cfg_in.groups.product_group.processing_type + + # product version + product_version_float = cfg_in.groups.product_group.product_version + product_version = f'{product_version_float:.1f}' + + # DEM description + dem_description = cfg_in.dem_description + + if not dem_description: + # If the DEM description is not provided, use DEM source + dem_description = os.path.basename(cfg_in.dem) # Manifests the field names, corresponding values from RTC workflow, and the description. # To extend this, add the lines with the format below: @@ -141,12 +209,16 @@ def populate_metadata_group(h5py_obj: h5py.File, # 'identification/relativeOrbitNumber': # [int(burst_in.burst_id[1:4]), 'Relative orbit number'], 'identification/trackNumber': - [int(burst_in.burst_id.split('_')[1]), 'Track number'], + [burst_in.burst_id.track_number, 'Track number'], + 'identification/burstID': + [str(burst_in.burst_id), 'Burst identification (burst ID)'], + 'identification/boundingPolygon': + [get_polygon_wkt(burst_in), + 'OGR compatible WKT representation of bounding polygon of the image'], 'identification/missionId': [burst_in.platform_id, 'Mission identifier'], - # NOTE maybe `SLC` has to be sth. like RTC? 'identification/productType': - ['SLC', 'Product type'], + ['RTC-S1', 'Product type'], # NOTE: in NISAR, the value has to be in UPPERCASE or lowercase? 'identification/lookDirection': ['Right', 'Look direction can be left or right'], @@ -168,7 +240,9 @@ def populate_metadata_group(h5py_obj: h5py.File, 'identification/diagnosticModeFlag': [False, 'Indicates if the radar mode is a diagnostic mode or not: True or False'], 'identification/processingType': - ['UNDEFINED', 'NOMINAL (or) URGENT (or) CUSTOM (or) UNDEFINED'], + [processing_type, 'NOMINAL (or) URGENT (or) CUSTOM (or) UNDEFINED'], + 'identification/productVersion': + [product_version, 'Product version'], # 'identification/frameNumber': # TBD # 'identification/productVersion': # Defined by RTC SAS # 'identification/plannedDatatakeId': @@ -201,6 +275,10 @@ def populate_metadata_group(h5py_obj: h5py.File, 'Radiometric terrain correction (RTC) algorithm'], 'RTC/metadata/processingInformation/algorithms/ISCEVersion': [isce3.__version__, 'ISCE version used for processing'], + 'RTC/metadata/processingInformation/algorithms/RTCVersion': + [str(__version__), 'RTC-S1 SAS version used for processing'], + 'RTC/metadata/processingInformation/algorithms/S1ReaderVersion': + [release_version, 'S1-Reader version used for processing'], 'RTC/metadata/processingInformation/inputs/l1SlcGranules': [l1_slc_granules, 'List of input L1 RSLC products used'], @@ -211,8 +289,8 @@ def populate_metadata_group(h5py_obj: h5py.File, 'List of input calibration files used'], 'RTC/metadata/processingInformation/inputs/configFiles': [cfg_in.run_config_path, 'List of input config files used'], - 'RTC/metadata/processingInformation/inputs/demFiles': - [dem_files, 'List of input dem files used'] + 'RTC/metadata/processingInformation/inputs/demSource': + [dem_description, 'DEM source description'] } for fieldname, data in dict_field_and_data.items(): path_dataset_in_h5 = os.path.join(root_path, fieldname) diff --git a/src/rtc/helpers.py b/src/rtc/helpers.py index 5b5c7e44..5c7c041e 100644 --- a/src/rtc/helpers.py +++ b/src/rtc/helpers.py @@ -1,10 +1,15 @@ '''collection of useful functions used across workflows''' import os - -import isce3 +import sqlite3 import logging +import journal +import numpy as np +from pyproj.transformer import Transformer from osgeo import gdal +from shapely import geometry + +import isce3 import rtc @@ -150,3 +155,164 @@ def check_dem(dem_path: str): err_str = f'DEM epsg of {epsg} out of bounds' logger.error(err_str) raise ValueError(err_str) + + + +def bbox_to_utm(bbox, *, epsg_src, epsg_dst): + """Convert bounding box coordinates to UTM. + + Parameters + ---------- + bbox : tuple + Tuple containing the lon/lat bounding box coordinates + (left, bottom, right, top) in degrees + epsg_src : int + EPSG code identifying input bbox coordinate system + epsg_dst : int + EPSG code identifying output coordinate system + + Returns + ------- + tuple + Tuple containing the bounding box coordinates in UTM (meters) + (left, bottom, right, top) + """ + xmin, ymin, xmax, ymax = bbox + xys = _convert_to_utm([(xmin, ymin), (xmax, ymax)], epsg_src, epsg_dst) + return (*xys[0], *xys[1]) + + +def polygon_to_utm(poly, *, epsg_src, epsg_dst): + """Convert a shapely.Polygon's coordinates to UTM. + + Parameters + ---------- + poly: shapely.geometry.Polygon + Polygon object + epsg : int + EPSG code identifying output projection system + + Returns + ------- + tuple + Tuple containing the bounding box coordinates in UTM (meters) + (left, bottom, right, top) + """ + coords = np.array(poly.exterior.coords) + xys = _convert_to_utm(coords, epsg_src, epsg_dst) + return geometry.Polygon(xys) + + +def _convert_to_utm(points_xy, epsg_src, epsg_dst): + """Convert a list of points to a specified UTM coordinate system. + + If epsg_src is 4326 (lat/lon), assumes points_xy are in degrees. + """ + if epsg_dst == epsg_src: + return points_xy + + t = Transformer.from_crs(epsg_src, epsg_dst, always_xy=True) + xs, ys = np.array(points_xy).T + xt, yt = t.transform(xs, ys) + return list(zip(xt, yt)) + + +def burst_bbox_from_db(burst_id, burst_db_file=None, burst_db_conn=None): + """Find the bounding box of a burst in the database. + + Parameters + ---------- + burst_id : str + JPL burst ID + burst_db_file : str + Location of burst database sqlite file, by default None + burst_db_conn : sqlite3.Connection + Connection object to burst database (If already connected) + Alternative to providing burst_db_file, will be faster + for multiply queries. + + Returns + ------- + epsg : int + EPSG code(s) of burst bounding box(es) + bbox : tuple[float] + Bounding box of burst in EPSG coordinates. Bounding box given as + tuple(xmin, ymin, xmax, ymax) + + Raises + ------ + ValueError + If burst_id is not found in burst database + """ + # example burst db: + # /home/staniewi/dev/burst_map_IW_000001_375887.OPERA-JPL.sqlite3 + if burst_db_conn is None: + burst_db_conn = sqlite3.connect(burst_db_file) + burst_db_conn.row_factory = sqlite3.Row # return rows as dicts + + query = "SELECT epsg, xmin, ymin, xmax, ymax FROM burst_id_map WHERE burst_id_jpl = ?" + cur = burst_db_conn.execute(query, (burst_id,)) + result = cur.fetchone() + + if not result: + raise ValueError(f"Failed to find {burst_id} in {burst_db_file}") + + epsg = result["epsg"] + bbox = (result["xmin"], result["ymin"], result["xmax"], result["ymax"]) + + return epsg, bbox + + +def burst_bboxes_from_db(burst_ids, burst_db_file=None, burst_db_conn=None): + """Find the bounding box of bursts in the database. + + Parameters + ---------- + burst_id : list[str] + list of JPL burst IDs. + burst_db_file : str + Location of burst database sqlite file, by default None + burst_db_conn : sqlite3.Connection + Connection object to burst database (If already connected) + Alternative to providing burst_db_file, will be faster + for multiply queries. + + Returns + ------- + bboxes : dict + Burst bounding boxes as a dict with burst IDs as key and tuples of + EPSG and bounding boxes (tuple[float]) as values. Bounding box given as + tuple(xmin, ymin, xmax, ymax) + + Raises + ------ + ValueError + If no burst_ids are found in burst database + """ + # example burst db: + # /home/staniewi/dev/burst_map_IW_000001_375887.OPERA-JPL.sqlite3 + if burst_db_conn is None: + burst_db_conn = sqlite3.connect(burst_db_file) + burst_db_conn.row_factory = sqlite3.Row # return rows as dicts + + # concatenate '?, ' with for each burst ID for IN query + qs_in_query = ', '.join('?' for _ in burst_ids) + query = f"SELECT * FROM burst_id_map WHERE burst_id_jpl IN ({qs_in_query})" + cur = burst_db_conn.execute(query, burst_ids) + results = cur.fetchall() + + if not results: + raise ValueError(f"Failed to find {burst_ids} in {burst_db_file}") + + n_results = len(results) + epsgs = [[]] * n_results + bboxes = [[]] * n_results + burst_ids = [[]] * n_results + for i_result, result in enumerate(results): + epsgs[i_result] = result["epsg"] + bboxes[i_result] = (result["xmin"], result["ymin"], + result["xmax"], result["ymax"]) + burst_ids[i_result] = result["burst_id_jpl"] + + # TODO add warning if not all burst bounding boxes found + return dict(zip(burst_ids, zip(epsgs, bboxes))) diff --git a/src/rtc/mosaic_geobursts.py b/src/rtc/mosaic_geobursts.py index 1fe9eaef..3f32d591 100644 --- a/src/rtc/mosaic_geobursts.py +++ b/src/rtc/mosaic_geobursts.py @@ -103,22 +103,22 @@ def check_reprojection(geogrid_mosaic, return flag_requires_reprojection -def weighted_mosaic(list_rtc_images, list_nlooks, geo_filename, - geogrid_in=None, verbose = True): +def compute_weighted_mosaic_array(list_rtc_images, list_nlooks, + geogrid_in=None, verbose = True): ''' - Mosaic the snapped S1 geobursts + Mosaic S-1 geobursts and return the mosaic as dictionary paremeters: ----------- list_rtc: list - list of the path to the rtc geobursts + List of the path to the rtc geobursts list_nlooks: list - list of the nlooks raster that corresponds to list_rtc - geo_filename: str - Path to the output mosaic + List of the nlooks raster that corresponds to list_rtc geogrid_in: isce3.product.GeoGridParameters, default: None - geogrid information to determine the output mosaic's shape and projection + Geogrid information to determine the output mosaic's shape and projection The geogrid of the output mosaic will automatically determined when it is None - + Returns: + mosaic_dict: dict + Mosaic dictionary ''' num_raster = len(list_rtc_images) @@ -238,15 +238,69 @@ def weighted_mosaic(list_rtc_images, list_nlooks, geo_filename, raster_rtc = None raster_nlooks = None - # Retreive the datatype information from the first input image + for i_band in range(num_bands): + valid_ind = np.where(arr_denominator > 0) + arr_numerator[i_band][valid_ind] = \ + arr_numerator[i_band][valid_ind] / arr_denominator[valid_ind] + + invalid_ind = np.where(arr_denominator == 0) + arr_numerator[i_band][invalid_ind] = np.nan + + mosaic_dict = { + 'mosaic_array': arr_numerator, + 'length': dim_mosaic[0], + 'width': dim_mosaic[1], + 'num_bands': num_bands, + 'wkt_projection': wkt_projection, + 'xmin_mosaic': xmin_mosaic, + 'ymax_mosaic': ymax_mosaic, + 'posting_x': posting_x, + 'posting_y': posting_y + } + return mosaic_dict + + + + +def compute_weighted_mosaic_raster(list_rtc_images, list_nlooks, geo_filename, + geogrid_in=None, verbose = True): + ''' + Mosaic the snapped S1 geobursts + paremeters: + ----------- + list_rtc: list + List of the path to the rtc geobursts + list_nlooks: list + List of the nlooks raster that corresponds to list_rtc + geo_filename: str + Path to the output mosaic + geogrid_in: isce3.product.GeoGridParameters, default: None + Geogrid information to determine the output mosaic's shape and projection + The geogrid of the output mosaic will automatically determined when it is None + + ''' + mosaic_dict = compute_weighted_mosaic_array(list_rtc_images, list_nlooks, + geogrid_in=geogrid_in, verbose = verbose) + + arr_numerator = mosaic_dict['mosaic_array'] + length = mosaic_dict['length'] + width = mosaic_dict['width'] + num_bands = mosaic_dict['num_bands'] + wkt_projection = mosaic_dict['wkt_projection'] + xmin_mosaic = mosaic_dict['xmin_mosaic'] + ymax_mosaic = mosaic_dict['ymax_mosaic'] + posting_x = mosaic_dict['posting_x'] + posting_y = mosaic_dict['posting_y'] + + # Retrieve the datatype information from the first input image reference_raster = gdal.Open(list_rtc_images[0], gdal.GA_ReadOnly) datatype_mosaic = reference_raster.GetRasterBand(1).DataType reference_raster = None - # write out the array + # Write out the array drv_out = gdal.GetDriverByName('Gtiff') raster_out = drv_out.Create(geo_filename, - dim_mosaic[1], dim_mosaic[0], num_bands, + width, length, num_bands, datatype_mosaic) raster_out.SetGeoTransform((xmin_mosaic, posting_x, 0, ymax_mosaic, 0, posting_y)) @@ -254,11 +308,63 @@ def weighted_mosaic(list_rtc_images, list_nlooks, geo_filename, raster_out.SetProjection(wkt_projection) for i_band in range(num_bands): - valid_ind = np.where(arr_denominator > 0) - arr_numerator[i_band][valid_ind] = \ - arr_numerator[i_band][valid_ind] / arr_denominator[valid_ind] + raster_out.GetRasterBand(i_band+1).WriteArray(arr_numerator[i_band]) - invalid_ind = np.where(arr_denominator == 0) - arr_numerator[i_band][invalid_ind] = np.nan - raster_out.GetRasterBand(i_band+1).WriteArray(arr_numerator[i_band]) + +def compute_weighted_mosaic_raster_single_band(list_rtc_images, list_nlooks, + output_file_list, + geogrid_in=None, verbose = True): + ''' + Mosaic the snapped S1 geobursts + paremeters: + ----------- + list_rtc: list + List of the path to the rtc geobursts + list_nlooks: list + List of the nlooks raster that corresponds to list_rtc + output_file_list: list + Output file list + geogrid_in: isce3.product.GeoGridParameters, default: None + Geogrid information to determine the output mosaic's shape and projection + The geogrid of the output mosaic will automatically determined when it is None + + ''' + mosaic_dict = compute_weighted_mosaic_array(list_rtc_images, list_nlooks, + geogrid_in=geogrid_in, verbose = verbose) + + arr_numerator = mosaic_dict['mosaic_array'] + length = mosaic_dict['length'] + width = mosaic_dict['width'] + num_bands = mosaic_dict['num_bands'] + wkt_projection = mosaic_dict['wkt_projection'] + xmin_mosaic = mosaic_dict['xmin_mosaic'] + ymax_mosaic = mosaic_dict['ymax_mosaic'] + posting_x = mosaic_dict['posting_x'] + posting_y = mosaic_dict['posting_y'] + + if num_bands != len(output_file_list): + error_str = (f'ERROR number of output files ({len(output_file_list)})' + f' does not match with the number' + f' of input bursts` bands ({num_bands})') + raise ValueError(error_str) + + for i_band, output_file in enumerate(output_file_list): + + # Retrieve the datatype information from the first input image + reference_raster = gdal.Open(list_rtc_images[0], gdal.GA_ReadOnly) + datatype_mosaic = reference_raster.GetRasterBand(1).DataType + reference_raster = None + + # Write out the array + drv_out = gdal.GetDriverByName('Gtiff') + raster_out = drv_out.Create(output_file, + width, length, num_bands, + datatype_mosaic) + + raster_out.SetGeoTransform((xmin_mosaic, posting_x, 0, ymax_mosaic, 0, posting_y)) + + raster_out.SetProjection(wkt_projection) + + for i_band in range(num_bands): + raster_out.GetRasterBand(i_band+1).WriteArray(arr_numerator[i_band]) diff --git a/src/rtc/runconfig.py b/src/rtc/runconfig.py index d7f7b3cf..7b6aefca 100644 --- a/src/rtc/runconfig.py +++ b/src/rtc/runconfig.py @@ -13,7 +13,7 @@ from rtc import helpers from rtc.radar_grid import file_to_rdr_grid -from rtc.geogrid import generate_geogrids +from rtc.geogrid import generate_geogrids, generate_geogrids_from_db from rtc.wrap_namespace import wrap_namespace, unwrap_to_dict from s1reader.s1_burst_slc import Sentinel1BurstSlc from s1reader.s1_orbit import get_orbit_file_from_list @@ -129,11 +129,11 @@ def validate_group_dict(group_cfg: dict, workflow_name) -> None: helpers.check_file_path(dem_path) helpers.check_dem(dem_path) - # Check 'product_path_group' section of runconfig. + # Check 'product_group' section of runconfig. # Check that directories herein have writing permissions - product_path_group = group_cfg['product_path_group'] - helpers.check_write_dir(product_path_group['product_path']) - helpers.check_write_dir(product_path_group['scratch_path']) + product_group = group_cfg['product_group'] + helpers.check_write_dir(product_group['product_path']) + helpers.check_write_dir(product_group['scratch_path']) def runconfig_to_bursts(cfg: SimpleNamespace) -> list[Sentinel1BurstSlc]: @@ -193,7 +193,7 @@ def runconfig_to_bursts(cfg: SimpleNamespace) -> list[Sentinel1BurstSlc]: for burst in load_bursts(safe_file, orbit_path, i_subswath, pol, flag_apply_eap=False): # get burst ID - burst_id = burst.burst_id + burst_id = str(burst.burst_id) # is burst_id wanted? skip if not given in config if (not cfg.input_file_group.burst_id is None and @@ -334,9 +334,14 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str) -> RunConfig: # Load geogrids dem_file = groups_cfg['dynamic_ancillary_file_group']['dem_file'] - geogrid_all, geogrids = generate_geogrids(bursts, geocoding_dict, - dem_file) + burst_database_file = groups_cfg['static_ancillary_file_group']['burst_database_file'] + if burst_database_file is None: + geogrid_all, geogrids = generate_geogrids(bursts, geocoding_dict, dem_file) + else: + geogrid_all, geogrids = generate_geogrids_from_db(bursts, geocoding_dict, + dem_file, burst_database_file) + # Empty reference dict for base runconfig class constructor empty_ref_dict = {} @@ -355,6 +360,10 @@ def burst_id(self) -> list[str]: def dem(self) -> str: return self.groups.dynamic_ancillary_file_group.dem_file + @property + def dem_description(self) -> str: + return self.groups.dynamic_ancillary_file_group.dem_description + @property def is_reference(self) -> bool: return self.groups.input_file_group.reference_burst.is_reference @@ -369,7 +378,7 @@ def polarization(self) -> list[str]: @property def product_path(self): - return self.groups.product_path_group.product_path + return self.groups.product_group.product_path @property def reference_path(self) -> str: @@ -397,11 +406,11 @@ def safe_files(self) -> list[str]: @property def product_id(self): - return self.groups.product_path_group.product_id + return self.groups.product_group.product_id @property def scratch_path(self): - return self.groups.product_path_group.scratch_path + return self.groups.product_group.scratch_path @property def gpu_enabled(self): @@ -425,7 +434,7 @@ def as_dict(self): date_str = lambda b : b.sensing_start.date().strftime('%Y%m%d') # create an unique burst key - burst_as_key = lambda b : '_'.join([b.burst_id, + burst_as_key = lambda b : '_'.join([str(b.burst_id), date_str(b), b.polarization]) diff --git a/src/rtc/schemas/rtc_s1.yaml b/src/rtc/schemas/rtc_s1.yaml index 36415e26..02c54c98 100644 --- a/src/rtc/schemas/rtc_s1.yaml +++ b/src/rtc/schemas/rtc_s1.yaml @@ -3,50 +3,73 @@ runconfig: groups: pge_name_group: - pge_name: enum('RTC_S1_PGE') + pge_name: enum('RTC_S1_PGE') input_file_group: - # Required. List of SAFE files (min=1) - safe_file_path: list(str(), min=1) - # Required. List of orbit (EOF) files - orbit_file_path: list(str(), min=1) - # Required. The unique burst ID to process - burst_id: list(str(), min=1, required=False) + # Required. List of SAFE files (min=1) + safe_file_path: list(str(), min=1) + # Required. List of orbit (EOF) files + orbit_file_path: list(str(), min=1) + # Optional. Burst ID to process (empty for all bursts) + burst_id: list(str(), min=1, required=False) dynamic_ancillary_file_group: - # Digital Elevation Model. - dem_file: str(required=False) - - product_path_group: - # Directory where PGE will place results - product_path: str() - # Directory where SAS can write temporary data - scratch_path: str() - - # If option `mosaic_bursts` is not set, output files are saved to: - # {output_dir}/{burst_id}/{product_id}{suffix}.{ext} - # If option `mosaic_bursts` is set, output files are saved to: - # {output_dir}/{product_id}{suffix}.{ext} - # If the field `product_id`` is left empty, the prefix "rtc_product" - # will be used instead. - # `suffix` is only used when there are multiple output files. - # `ext` is determined by geocoding_options.output_format. - output_dir: str() - product_id: str() - mosaic_bursts: bool(required=False) - - # Format of output file - output_format: enum('HDF5', 'NETCDF', 'ENVI', 'GTiff', 'COG', required=False) + # Digital Elevation Model. + dem_file: str(required=False) + + # Digital elevation model description + dem_description: str(required=False) + + static_ancillary_file_group: + + # burst database sqlite file + burst_database_file: str(required=False) primary_executable: - product_type: enum('RTC_S1') + product_type: enum('RTC_S1') + + product_group: + + processing_type: enum('NOMINAL', 'URGENT', 'CUSTOM', 'UNDEFINED', required=False) + + product_version: num(required=False) + + # Directory where PGE will place results + product_path: str() + # Directory where SAS can write temporary data + scratch_path: str() + + # If option `save_bursts` is set, output bursts are saved to: + # {output_dir}/{burst_id}/{product_id}_v{product_version}{suffix}.{ext} + # If option `save_mosaics` is set, output mosaics are saved to: + # {output_dir}/{product_id}_v{product_version}{suffix}.{ext} + # If the field `product_id`` is left empty, the prefix "rtc_product" + # will be used instead. + # `suffix` is only used when there are multiple output files. + # `ext` is determined by geocoding_options.output_imagery_format. + output_dir: str() + product_id: str() + + # RTC-S1 imagery + save_bursts: bool(required=False) + save_mosaics: bool(required=False) + output_imagery_format: enum('HDF5', 'NETCDF', 'ENVI', 'GTiff', 'COG', required=False) + output_imagery_compression: str(required=False) + output_imagery_nbits: int(min=1, required=False) + + # Optional. Save secondary layers (e.g., inc. angle) within + # the HDF5 file + save_secondary_layers_as_hdf5: bool(required=False) + + # Save RTC-S1 metadata in the HDF5 format + # Optional for `output_imagery_format` equal to 'ENVI', 'GTiff', or + # 'COG', and enabled by default for `output_imagery_format` equal + # to 'HDF5' or 'NETCDF' or `save_secondary_layers_as_hdf5` is True + save_metadata: bool(required=False) # This section includes parameters to tweak the workflow processing: include('processing_options', required=False) - # Worker options (e.g. enable/disable GPU processing, select GPU device ID) - worker: include('worker_options', required=False) - --- geo2rdr_options: @@ -55,16 +78,31 @@ geo2rdr_options: # Maximum number of iterations numiter: int(min=1, required=False) +rdr2geo_options: + # Convergence threshold for rdr2geo algorithm + threshold: num(min=0, required=False) + # Maximum number of iterations + numiter: int(min=1, required=False) + # Group of processing options processing_options: - # Polarization to process. 3 modes below correspond to VV, VH, VV+VH + + # Check if ancillary inputs cover entirely the output product + check_ancillary_inputs_coverage: bool(required=False) + + # Polarization channels to process. 3 modes below correspond to VV, VH, VV+VH polarization: enum('co-pol', 'cross-pol', 'dual-pol', required=False) # Options to run geo2rdr geo2rdr: include('geo2rdr_options', required=False) + + # Options to run rdr2geo (for running topo when calculating layover shadow mask) + rdr2geo: include('rdr2geo_options', required=False) + # Range split-spectrum options range_split_spectrum: include('range_split_spectrum_options', required=False) + # DEM interpolation method dem_interpolation_method: enum('sinc', 'bilinear', 'bicubic', 'nearest', 'biquintic', required=False) # Apply absolute radiometric correction @@ -76,6 +114,12 @@ processing_options: # Apply RTC apply_rtc: bool(required=False) + # Apply bistatic delay correction + apply_bistatic_delay_correction: bool(required=False) + + # Apply dry tropospheric delay correction + apply_dry_tropospheric_delay_correction: bool(required=False) + # Radiometric Terrain Correction (RTC) rtc: include('rtc_options', required=False) @@ -136,6 +180,9 @@ geocoding_options: # Save the interpolated DEM used to generate the RTC product save_dem: bool(required=False) + # Save layover shadow mask + save_layover_shadow_mask: bool(required=False) + # OPTIONAL - Absolute radiometric correction abs_rad_cal: num(required=False)