Merge pull request #32 from EarthObservation/save_vrt

Save vrt

Author: zkokalj

README.md

@@ -1,6 +1,6 @@
 [![PyPI](https://img.shields.io/pypi/v/RVT_py?style=flat-square)](https://pypi.org/project/rvt-py/)
-[![Anaconda-Server Badge](https://anaconda.org/rvtpy/rvt_py/badges/version.svg)](https://anaconda.org/rvtpy/rvt_py)
-[![Anaconda-Server Badge](https://anaconda.org/rvtpy/rvt_py/badges/latest_release_date.svg)](https://anaconda.org/rvtpy/rvt_py)
+[![Anaconda-Server Badge](https://anaconda.org/zmigyyy/rvt_py/badges/version.svg)](https://anaconda.org/zmigyyy/rvt_py)
+[![Anaconda-Server Badge](https://anaconda.org/zmigyyy/rvt_py/badges/latest_release_date.svg)](https://anaconda.org/zmigyyy/rvt_py)
 
 # Relief Visualization Toolbox Python library
 
@@ -21,7 +21,7 @@ Methods currently implemented are:
 * sky-view factor (as developed by our team),
 * anisotropic sky-view factor,
 * positive and negative openness,
-* local dominance, 
+* local dominance,
 * multi-scale topographic position.
 
 ## RVT for Python
@@ -43,10 +43,11 @@ When using the tools, please cite:
 
 ## Installation
 
-The RVT Python package can be installed using Conda or PyPI, and can be used in Python scripts, Jupyter Notebooks, QGIS and ArcGIS Pro.
+The RVT Python package can be installed using Conda or PyPI, and can be used in Python scripts, Jupyter Notebooks and ArcGIS Pro.
 
-RVT can also be installed as [a set of custom raster functions for ArcGIS](https://rvt-py.readthedocs.io/en/latest/install_arcgis.html "ArcGIS installation") and [a plugin for QGIS](https://rvt-py.readthedocs.io/en/latest/install_qgis.html "QGIS installation").
+RVT can also be installed as [a set of custom raster functions for ArcGIS](https://rvt-py.readthedocs.io/en/latest/install_arcgis.html "ArcGIS installation"), and [a plugin for QGIS](https://rvt-py.readthedocs.io/en/latest/install_qgis.html "QGIS installation").
 
+You can also clone the repository.
 
 ### Conda
 
@@ -70,25 +71,14 @@ Then open Command Prompt (Windows) or Terminal (MacOS) and run:
 
 ### Requirements
 
-*   Python 3.8+
+Required libraries (specified versions have been tested, other versions may also work):
 
-> Code was tested on python version 3.11.11
+*   numpy 1.19.2
+*   scipy 1.5.2
+*   gdal 3.0.2
+*   rasterio 1.2.6
 
-Dependencies for core functionality:
-
-*   gdal
-*   matplotlib
-*   numpy
-*   scipy
-
-Dependencies for extended use (running scripts in examples folder):
-
-*   pandas 
-*   geopandas
-*   jupyter
-*   rasterio
-
-> We recommend using `Anaconda` virtual environments and `conda-forge` channel for installing the required packages (this works best with ``gdal``).
+We recommend using Python 3.6 or higher and a Conda environment (this works best with ``gdal``).
 
 ## Documentation
 Documentation of the package and its use is available at [Relief Visualization Toolbox in Python documentation](https://rvt-py.readthedocs.io/).
@@ -103,7 +93,7 @@ Development of RVT Python scripts was part financed by the Slovenian Research Ag
 This project is licensed under the terms of the [Apache License](LICENSE).
 
 ## About
-RVT Python library by Žiga Kokalj, Žiga Maroh, Krištof Oštir, Klemen Zakšek and Nejc Čož, 2024.
+RVT Python library by Žiga Kokalj, Žiga Maroh, Krištof Oštir, Klemen Zakšek and Nejc Čož, 2022.
 
 It is developed in collaboration between ZRC SAZU and University of Ljubljana. 
 

examples/tiled_processing/_test_tiled.py

@@ -3,19 +3,18 @@
 if __name__ == "__main__":
     # List of INPUT DATASETS (have to be in this format)
     list_tifs = [
-        r"c:\Users\ncoz\GitHub\erc_potencial\test_data\test_small.tif",
-        # r"c:\test_data\testmx\2022MxCP_C_dem_05m.tif",
-        # r"c:\test_data\testmx_s\2022MxCP_C_dem_05m_s.tif"
+        r"c:\path-to-file\test-file.tif",
     ]
 
     # Select visualizations
     # 'slp', 'hs', 'slrm', 'svf', 'opns', 'neg_opns', 'ld', 'mstp'
     vis_types = [
-        'slrm'
+        "slrm"
     ]
     # "vat_combined" / "e2MSTP" / "e4MSTP" / "crim" / "rrim"
     blend_types = [
-        "vat_combined"
+        # "vat_combined"
     ]
 
-    run_main(list_tifs, vis_types, blend_types, save_float=False)
+    # save_vrt (if True saves VRT, if False saves as TIF mosaic), only used when tiling
+    run_main(list_tifs, vis_types, blend_types, save_float=False, save_vrt=True)

examples/tiled_processing/rvt_tiled.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "336b1fb4-8369-43b3-9b99-7915bad70954",
    "metadata": {},
    "outputs": [],
@@ -22,14 +22,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "cc536aa8-15f8-4e28-bbab-4c1c5c8f52a3",
    "metadata": {},
    "outputs": [],
    "source": [
     "# List of INPUT DATASETS (have to be in this format)\n",
     "list_tifs = [\n",
-    "    r\"c:\\Users\\ncoz\\GitHub\\erc_potencial\\test_data\\test_small.tif\"\n",
+    "    r\"c:\\test_data\\2022MxCP_1GB\\2022MxCP_E_dem_05m.tif\"\n",
     "]"
    ]
   },
@@ -38,7 +38,13 @@
    "id": "caac01bb-d3f1-45d9-b041-ca200f15f44d",
    "metadata": {},
    "source": [
-    "The defaults are save in 8-bit format (256 values, RGB). To also save pixel values as float, change the following parameter to \"True\":"
+    "**save_float**\n",
+    "* The defaults are saved in 8-bit format (256 values, RGB)\n",
+    "* To also save pixel values as float, change the `save_float` parameter to `True`\n",
+    "\n",
+    "**save_vrt** \n",
+    "* if `True` the results are left as tiles and a VRT file is created\n",
+    "* if `False` results are merged into TIF mosaic and indvidual tiles are deleted"
    ]
   },
   {
@@ -48,7 +54,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "save_float = False"
+    "save_float = False\n",
+    "save_vrt = True"
    ]
   },
   {
@@ -83,14 +90,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "id": "01d02ac4-16cc-45a9-88d4-ecea209e9a7f",
    "metadata": {},
    "outputs": [],
    "source": [
     "# Select visualizations\n",
     "vis_types = [\n",
-    "    \"slrm\"\n",
+    "    \n",
     "]"
    ]
   },
@@ -122,7 +129,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "id": "374da272-c812-4dbd-98c4-32e3839c6d60",
    "metadata": {},
    "outputs": [],
@@ -145,26 +152,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "id": "0e55d218-4923-440d-988b-75824aceeb72",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Start --- c:/Users/ncoz/GitHub/erc_potencial/test_data/test_small.tif\n",
-      "Processing as single tile...\n",
-      "Done with computing blends in 0 min.\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "run_main(\n",
     "    list_tifs=list_tifs,\n",
     "    vis_types=vis_types,\n",
     "    blend_types=blend_types,\n",
-    "    save_float=save_float\n",
+    "    save_float=save_float,\n",
+    "    save_vrt=save_vrt\n",
     ")"
    ]
   }

examples/tiled_processing/rvt_tiled_GUI.ipynb

@@ -62,25 +62,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "4b90a3fd-1fe7-46d4-9cab-e96137dfdd48",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "910520b1f2ea480bbf18ec0a19469004",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "VBox(children=(HTML(value='<b>RVT - TILED PROCESSING FOR LARGE FILES</b>'), SelectFilesButton(description='Sel…"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "import tiled_GUI"
    ]

examples/tiled_processing/tiled_GUI.py

@@ -126,6 +126,7 @@ def on_button_clicked(b):
         # print(str(save_float_checkbox.value))
 
         run_main(list_tifs, vis_types, blend_types, save_float)
+        # todo: add save_vrt=True to arguments and create a checkbox for this
 
     button_run_adaf.disabled = False
 

examples/tiled_processing/tiled_multiprocess.py

@@ -13,9 +13,11 @@
 
 import numpy as np
 import rasterio
+from geopandas import GeoDataFrame
 from osgeo import gdal
 from rasterio.merge import merge
 from rasterio.windows import from_bounds
+from shapely.geometry import box
 
 import grid_tools as gt
 import rvt.blend
@@ -26,7 +28,7 @@
 gdal.UseExceptions()
 
 
-def run_main(list_tifs, vis_types, blend_types, save_float=False):
+def run_main(list_tifs, vis_types, blend_types, save_float=False, save_vrt=False):
     for in_file in list_tifs:
         in_file = Path(in_file)
 
@@ -44,6 +46,7 @@ def run_main(list_tifs, vis_types, blend_types, save_float=False):
                 tile_size = 2000
             else:
                 tile_size = None
+            raster_bounds = src.bounds
 
         if tile_size:
             # # (3) To filter we need polygon covering valid data
@@ -62,6 +65,10 @@ def run_main(list_tifs, vis_types, blend_types, save_float=False):
             #     save_gpkg=False
             # )
 
+            # Crop outer tiles to the extents of raster (IMPORTANT so output size is same as input when tiling!)
+            raster_bbox = box(*raster_bounds)
+            tiles_extents = GeoDataFrame(geometry=tiles_extents.geometry.intersection(raster_bbox), crs=tiles_extents.crs)
+
             # Add extents column: Extents are (L, B, R, T).
             tiles_extents["extents"] = tiles_extents.bounds.apply(lambda x: (x.minx, x.miny, x.maxx, x.maxy), axis=1)
 
@@ -76,11 +83,12 @@ def run_main(list_tifs, vis_types, blend_types, save_float=False):
             blend_types=blend_types,
             input_vrt_path=in_file,
             tiles_list=tiles_list,
-            save_float=save_float
+            save_float=save_float,
+            save_vrt=save_vrt
         )
 
 
-def tiled_blending(vis_types, blend_types, input_vrt_path, tiles_list, save_float):
+def tiled_blending(vis_types, blend_types, input_vrt_path, tiles_list, save_float, save_vrt):
     t0 = time.time()
 
     # Prepare paths
@@ -146,15 +154,17 @@ def tiled_blending(vis_types, blend_types, input_vrt_path, tiles_list, save_floa
                 sf = False
             mosaic_path = build_vrt(result_paths, vrt_path, save_float=sf)
 
-            # Create mosaic and crop to the original extents
-            vrt_to_mosaic(
-                vrt_path=mosaic_path,  # ll_path / "test.tif",
-                output_mosaic_path=ll_path / f"{result}.tif",
-            )
+            # Build a mosaic from VRT if selected #todo: before app starts give warning about size if VRT is not ON
+            if not save_vrt:
+                # Create mosaic and crop to the original extents
+                vrt_to_mosaic(
+                    vrt_path=mosaic_path,  # ll_path / "test.tif",
+                    output_mosaic_path=ll_path / f"{result}.tif",
+                )
 
-            # Delete temp files
-            vrt_path.unlink()
-            shutil.rmtree(result_parent)
+                # Delete VRT files
+                vrt_path.unlink()
+                shutil.rmtree(result_parent)
 
             # print(results)
 
@@ -1284,6 +1294,9 @@ def get_required_arrays(vis_types, blend_types):
         req_arrays["opns_1"] = True
         req_arrays["neg_opns_1"] = True
 
+    if "new_blend" in blend_types:
+        req_arrays["slrm_1"] = True
+
     return req_arrays
 
 
@@ -1490,6 +1503,7 @@ def get_raster_vrt(vrt_path, extents, buffer):
         vrt_nodata = vrt.nodata
         vrt_transform = vrt.transform
         vrt_crs = vrt.crs
+        nodata_val = vrt.nodata
 
         # ADD BUFFER TO EXTENTS (LBRT) - transform pixels to meters!
         buffer_m = buffer * vrt_res[0]
@@ -1508,8 +1522,13 @@ def get_raster_vrt(vrt_path, extents, buffer):
         # boundless - if window falls out of bounds, read it and fill with NaNs
         win_array = vrt.read(window=buff_window, boundless=True)
 
+        # Deal with nodata
+        if nodata_val is not None and not np.isnan(nodata_val):
+            # Ensure array is float to support np.nan
+            win_array = win_array.astype("float32", copy=False)
+            win_array[win_array == nodata_val] = np.nan
+
         # Save transform object of both extents (original and buffered)
-        buff_transform = vrt.window_transform(buff_window)
         orig_transform = vrt.window_transform(orig_window)
 
     # For raster with only one band, remove first axis from the array (RVT requirement)
@@ -1519,7 +1538,7 @@ def get_raster_vrt(vrt_path, extents, buffer):
     # Prepare output metadata profile
     out_profile = {
         'driver': 'GTiff',
-        'nodata': None,
+        'nodata': np.nan,
         'width':  win_array.shape[1] - 2 * buffer,
         'height':  win_array.shape[0] - 2 * buffer,
         'count':  1,
@@ -1532,7 +1551,6 @@ def get_raster_vrt(vrt_path, extents, buffer):
         "array": win_array,
         "resolution": vrt_res,
         "no_data": vrt_nodata,
-        "buff_transform": buff_transform,
         "orig_transform": orig_transform,
         "crs": vrt_crs,
         "profile": out_profile