# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernel_info:
#     name: srspy
#   kernelspec:
#     display_name: rscube
#     language: python
#     name: rscube
# ---

# The goal of this notebook is to create a datacube across different rasters of remotely sensed data with different coordinate reference systems and resolutions. Specifically, the new set of rasters will all have the same size and each pixel across the rasters will all correspond to the same geographic area.
#
# All the data over the region of interest should be placed in `data/tifs`. In this notebook, we will select a subset the area using one of the tifs and reproject all the other rasters into the same subset with identical `crs` and `geotransform`.
#
# This notebook is *optional* and can be used to apply the methods found here to different study areas. We include the reprojected datacube for subsequent analyses.

import rasterio
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
from tqdm import tqdm
from rscube import get_cropped_profile, reproject_arr_to_match_profile

# # Data Paths

data_path_name = 'data/tifs'
data_path = Path(data_path_name)
data_path.exists()

tif_paths = sorted(list(data_path.glob('*.tif')))
tif_paths

dest_dir_name = 'data/tif_datacube'
dest_dir = Path(dest_dir_name)
dest_dir.mkdir(exist_ok=True, parents=True)

# # Inspect the Data

with rasterio.open(str(tif_paths[0])) as ds:
    hh = ds.read(1)
    profile = ds.profile.copy()
profile

plt.imshow(hh, vmin=0, vmax=.5)

# We are going to reproject the data into a subset.

# +
sy = np.s_[1000: 2500]
sx = np.s_[3000: 5000]

plt.imshow(hh[sy, sx], vmin=0, vmax=.5)
# -

cropped_profile = get_cropped_profile(profile, sx, sy)
cropped_profile


# # Reproject All

# We write a function to reproject one raster and then apply to each raster in our folder. We use the function `reproject_arr_to_match_profile` which takes the metadata stored in the `rasterio` profile and makes sure the other rasters match that geographic metadata. The rest of the function is io and making sure the new rasters are saved to the desired location.

def reproject_one(path):
    with rasterio.open(path) as ds:
        arr = ds.read()
        src_profile = ds.profile
    
    arr_r, profile_r = reproject_arr_to_match_profile(arr, src_profile, cropped_profile)
    
    dest_path = dest_dir/path.name
    with rasterio.open(dest_path, 'w', **profile_r) as ds:
        ds.write(arr_r)

    return dest_path


# +
# We have this to make sure the above function works as expected
reproject_one(tif_paths[0])

# This is equivalent to
# [reproject_path(path) for path in tif_paths]
list(map(reproject_one, tqdm(tif_paths)))
# -



StarcoderdataJupyterScriptsDedupFiltered