<jupyter_start><jupyter_text>In this notebook, we demonstrate how to:1. reproject the ASF radiometric and terrain corrected data into a single coordinate reference system.2. crop the data to a smaller area3. filter the time series according to month if relevant<jupyter_code>import rasterio
from pathlib import Path
from tqdm import tqdm
import lxml.etree as etree
import datetime
import numpy as np
import matplotlib.pyplot as plt
from rscube import (get_cropped_profile,
                    reproject_arr_to_match_profile,
                    reproject_arr_to_match_profile)<jupyter_output><empty_output><jupyter_text>Setting Data PathsYou should have all the data downloaded using the example from the previous notebook in the directory:```./data/asf_data/```<jupyter_code>DATA_DIR_NAME = f'data/asf_data/'
DATA_DIR = Path(DATA_DIR_NAME)
DATA_DIR.exists()
hh_paths = sorted(list(DATA_DIR.glob('*/*HH.tif')))
hv_paths = sorted(list(DATA_DIR.glob('*/*HV.tif')))
dem_paths = sorted(list(DATA_DIR.glob('*/*dem*.tif')))
metadata_paths = sorted(list(DATA_DIR.glob('*/*iso.xml')))
metadata_paths, dem_paths<jupyter_output><empty_output><jupyter_text>Inspect Profiles<jupyter_code>def get_profile(path):
    with rasterio.open(path) as ds:
        p = ds.profile
    return p
profiles = [get_profile(p) for p in hh_paths]
profiles<jupyter_output><empty_output><jupyter_text>Get DatesThe ASF data is distributed with an xml file and we have determined which item contains the date and we extract that. This will work with ASF RTC ALOS-1 images only. Each sensor will have data stored in a different file and in a different format.<jupyter_code>def get_alos_date(metadata_xml_path):
    tree = etree.parse(open(metadata_xml_path))
    root = tree.getroot()

    dataAcquisition_elements = root.xpath('//gml:beginPosition', namespaces=root.nsmap)
    assert(len(dataAcquisition_elements) == 1)
    element = dataAcquisition_elements[0].text
    date = datetime.date(int(element[:4]), int(element[5:7]), int(element[8:10]))
    return date
dates = list(map(get_alos_date, metadata_paths))
dates<jupyter_output><empty_output><jupyter_text>Filter by DateWe only want the spring and summer months (April - September).<jupyter_code># This means we only reserve months 4, ..., 9
MONTHS_TO_USE = list(range(4, 10))
MONTHS_TO_USE
indices = [k for k, date in enumerate(dates) if date.month in MONTHS_TO_USE]
indices
metadata_paths_filtered = [metadata_paths[k] for k in indices]
dem_paths_filtered = [dem_paths[k] for k in indices]
hv_paths_filtered = [hv_paths[k] for k in indices]
hh_paths_filtered = [hh_paths[k] for k in indices]
profiles = [get_profile(p) for p in hh_paths_filtered]
profiles<jupyter_output><empty_output><jupyter_text>Get Cropped SubsetWe are going to view one area and get a subset and reproject all the other images into that frame.<jupyter_code>with rasterio.open(hv_paths[1]) as ds:
    X = ds.read(1)
    mask = ~ds.read_masks(1).astype(bool)
    X[mask] = np.nan
sy = np.s_[250:]
sx = np.s_[:1200]

plt.figure(figsize=(10, 10))
plt.imshow(10* np.log10(X[sy, sx]))
plt.colorbar(label='$\gamma^0$')
with rasterio.open(hv_paths[0]) as ds:
    PROFILE = ds.profile

CROPPED_PROFILE = get_cropped_profile(PROFILE, sx, sy)
CROPPED_PROFILE['nodata'] = np.nan
CROPPED_PROFILE
REFERENCE_PROFILE = CROPPED_PROFILE.copy()
REFERENCE_PROFILE['nodata'] = np.nan
def read_one(tif_path):
    with rasterio.open(tif_path) as ds:
        img = ds.read(1)
        profile = ds.profile
        mask = (img == 0)
        img[mask] = np.nan
    profile['nodata'] = np.nan
    return img, profile

def reproject_one(tif_path):
    img_src, profile_src = read_one(tif_path)
    
    img_r, profile_r = reproject_arr_to_match_profile(img_src, 
                                                      profile_src, 
                                                      REFERENCE_PROFILE, 
                                                      resampling='bilinear')
    # Clips image to between -30 and 0 db.
    img_r = np.clip(img_r, 1e-3, 1)
    
    return img_r, profile_r

def write_reprojected(img_r, original_tif_path, metadata_path, pol):
    date = get_alos_date(metadata_path)

    dest_dir = Path(f'{DATA_DIR}_reprojected/{pol}')
    dest_dir.mkdir(exist_ok=True, parents=True)
    dest_path = dest_dir/f'ALOS1_RTC_{pol}_{date.year}{date.month:02d}{date.day:02d}.tif'

    with rasterio.open(dest_path, 'w', **REFERENCE_PROFILE) as ds:
        ds.write(img_r.astype(REFERENCE_PROFILE['dtype']))
    return dest_path
reprojected_data_hh = list(map(reproject_one, tqdm(hh_paths_filtered)))
reprojected_images_hh, _ = zip(*reprojected_data_hh)
reprojected_data_hv = list(map(reproject_one, tqdm(hv_paths_filtered)))
reprojected_images_hv, _ = zip(*reprojected_data_hv)
from functools import reduce
def find_mask(img_0, img_1):
    return np.isnan(img_0) | np.isnan(img_1)

ts_mask = list(reduce(find_mask, reprojected_images_hh + reprojected_images_hv))[0]
plt.imshow(ts_mask)
def apply_mask(img):
    img_ = img.copy()
    (img_[0,...])[ts_mask] = np.nan
    return img_
reprojected_images_hh = list(map(apply_mask, tqdm(reprojected_images_hh)))
reprojected_images_hv = list(map(apply_mask, tqdm(reprojected_images_hv)))
from itertools import starmap

N = len(metadata_paths_filtered)
list(starmap(write_reprojected, tqdm(zip(reprojected_images_hh, hh_paths_filtered, metadata_paths_filtered, ['hh'] * N), total=N)))
from itertools import starmap

N = len(metadata_paths_filtered)
list(starmap(write_reprojected, tqdm(zip(reprojected_images_hv, hv_paths_filtered, metadata_paths_filtered, ['hv'] * N), total=N)))<jupyter_output>100%|██████████| 4/4 [00:00<00:00, 14.55it/s]<jupyter_text>Cropping the DEMWe'll copy the dem from the first RTC image into this directory because this will be useful for filtering by slope.<jupyter_code>dest_dir = Path(f'{DATA_DIR}_reprojected')
dest_dir.exists(), dest_dir
with rasterio.open(dem_paths_filtered[0]) as ds:
    dem = ds.read(1)
    dem = dem[sy, sx]
    
dem_profile = REFERENCE_PROFILE.copy()
dem_profile['dtype'] = 'float32'
with rasterio.open(dest_dir/'dem.tif', 'w', **dem_profile) as ds:
    ds.write(dem.astype(np.float32), 1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11207869698435736, 'fraction_numerical': 0.012901144976616674, 'mean_word_length': 4.769302325581395, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 5}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '343976', 'license': 'Apache-2.0', 'n_tokens_mistral': 2300, 'n_tokens_neox': 2191, 'n_words': 488, 'path': 'notebooks/change_detection/ALOS1_Mangroves_Malaysia/1 - Rename and Reproject Time Series.ipynb', 'repo_name': 'simard-landscape-lab/rscube', 'size': '6201'}