<jupyter_start><jupyter_text>Subset Data Stack <NAME>; University of Alaska Fairbanks  This notebook crops a directory of tiffs to a subset area of interest using an interactive Matplotlib plot of an image in your data stack.  Important Note about JupyterHub  Your JupyterHub server will automatically shutdown when left idle for more than 1 hour. Your notebooks will not be lost but you will have to restart their kernels and re-run them from the beginning. You will not be able to seamlessly continue running a partially run notebook.    0. Importing Relevant Python Packages  In this notebook we will use the following scientific library:     GDAL is a software library for reading and writing raster and vector geospatial data formats. It includes a collection of programs tailored for geospatial data processing. Most modern GIS systems (such as ArcGIS or QGIS) use GDAL in the background. Import the necesssary libraries and modules:<jupyter_code>%%capture
import os
os.environ["PROJ_LIB"] = '/opt/conda/share/proj/'
import glob
import json # for loads
import shutil

from osgeo import gdal
import pyproj 
import numpy as np

from IPython.display import Markdown
from IPython.display import display

from matplotlib.widgets import RectangleSelector
import matplotlib.pyplot as plt 
plt.rcParams.update({'font.size': 12})

from asf_notebook import new_directory
from asf_notebook import path_exists
from asf_notebook import remove_nan_filled_tifs
from asf_notebook import AOI_Selector
from asf_notebook import select_parameter
from asf_notebook import get_polarity_from_path
from asf_notebook import date_from_product_name
from asf_notebook import jupytertheme_matplotlib_format<jupyter_output><empty_output><jupyter_text>Setup interactive Matplotlib plotting inside the notebook:<jupyter_code>%%capture
jupytertheme_matplotlib_format()
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>Write functions to gather and print individual tiff paths:<jupyter_code>def get_tiff_paths(paths):
    tiff_paths = glob.glob(paths)
    tiff_paths.sort()
    return tiff_paths

def print_tiff_paths(tiff_paths):
    print("Tiff paths:")
    for p in tiff_paths:
        print(f"{p}\n")<jupyter_output><empty_output><jupyter_text>Enter the path to the directory holding your tiffs:<jupyter_code># /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs
while True:
    print("Enter the absolute path to the directory holding your tiffs.")
    tiff_dir = input()
    paths = f"{tiff_dir}/*.tif*"
    if os.path.exists(tiff_dir):
        tiff_paths = get_tiff_paths(paths)
        if len(tiff_paths) < 1:
            print(f"{tiff_dir} exists but contains no tifs.")
            print("You will not be able to proceed until tifs are prepared.")
        break
    else:
        print(f"\n{tiff_dir} does not exist.")
        continue<jupyter_output>Enter the absolute path to the directory holding your tiffs.
/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs<jupyter_text>Determine the path to the analysis directory containing the tiff directory:<jupyter_code>analysis_dir = os.path.dirname(tiff_dir)
print(analysis_dir)<jupyter_output>/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer<jupyter_text>Determine the UTM zone for your images. This assumes you have already reprojected any tiffs with errant UTM zones to a single predominant UTM zone, using the Prepare_Data_Stack_Hyp3 notebook.<jupyter_code>info = (gdal.Info(tiff_paths[0], options = ['-json']))
info = json.dumps(info)
info = (json.loads(info))['coordinateSystem']['wkt']
utm = info.split('ID')[-1].split(',')[1][0:-2]
print(f"UTM Zone: {utm}")
tiff_paths = get_tiff_paths(paths)
print_tiff_paths(tiff_paths)<jupyter_output>Tiff paths:
/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180619T110424_DVP_RTC30_G_gpufem_3794_VH.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180619T110424_DVP_RTC30_G_gpufem_3794_VV.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180623T224352_DVP_RTC30_G_gpufem_B6E7_VH.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180623T224352_DVP_RTC30_G_gpufem_B6E7_VV.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180701T110424_DVP_RTC30_G_gpufem_B835_VH.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180701T110424_DVP_RTC30_G_gpufem_B835_VV.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180705T224418_DVP_RTC30_G_gpufem_0C30_VH.tif

/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180705T224[...]<jupyter_text>Create a string containing paths to one image for each area represented in the stack:<jupyter_code>to_merge = {}
for pth in tiff_paths:
    info = (gdal.Info(pth, options = ['-json']))
    info = json.dumps(info)
    info = (json.loads(info))['wgs84Extent']['coordinates']
    
    coords = [info[0][0], info[0][3]]
    for i in range(0, 2):
        for j in range(0, 2):
            coords[i][j] = round(coords[i][j])
    str_coords = f"{str(coords[0])}{str(coords[1])}"
    if str_coords not in to_merge:
        to_merge.update({str_coords: pth})
merge_paths = ""
for pth in to_merge:
    merge_paths = f"{merge_paths} {to_merge[pth]}"
print(merge_paths)<jupyter_output>/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180619T110424_DVP_RTC30_G_gpufem_3794_VH.tif /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180623T224352_DVP_RTC30_G_gpufem_B6E7_VH.tif<jupyter_text>Merge the images for display in the Area-Of-Interest selector:<jupyter_code>full_scene = f"{analysis_dir}/full_scene.tif"
if os.path.exists(full_scene):
    os.remove(full_scene)
gdal_command = f"gdal_merge.py -o {full_scene} {merge_paths}"
!{gdal_command}<jupyter_output>0...10...20...30...40...50...60...70...80...90...100 - done.<jupyter_text>Subset The Tiffs   Create a Virtual Raster Stack:<jupyter_code>image_file = f"{analysis_dir}/raster_stack.vrt"
!gdalbuildvrt -separate $image_file -overwrite $full_scene<jupyter_output>0...10...20...30...40...50...60...70...80...90...100 - done.<jupyter_text>Convert the VRT into an array:<jupyter_code>img = gdal.Open(image_file)
rasterstack = img.ReadAsArray()<jupyter_output><empty_output><jupyter_text>Print the number of bands, pixels, and lines:<jupyter_code>print(img.RasterCount) # Number of Bands
print(img.RasterXSize) # Number of Pixels
print(img.RasterYSize) # Number of Lines<jupyter_output>1
10770
12743<jupyter_text>Create an AOI selector from an image in your raster stack:<jupyter_code>fig_xsize = 7.5
fig_ysize = 7.5
aoi = AOI_Selector(rasterstack, fig_xsize, fig_ysize)<jupyter_output><empty_output><jupyter_text>Gather and define projection details:<jupyter_code>geotrans = img.GetGeoTransform()
projlatlon = pyproj.Proj('EPSG:4326') # WGS84
projimg = pyproj.Proj(f'EPSG:{utm}')<jupyter_output><empty_output><jupyter_text>Write a function to convert the pixel, line coordinates from the AOI selector into geographic coordinates in the stack's EPSG projection:<jupyter_code>def geolocation(x, y, geotrans,latlon=True):
    ref_x = geotrans[0]+x*geotrans[1]
    ref_y = geotrans[3]+y*geotrans[5]
    if latlon:
        ref_y, ref_x = pyproj.transform(projimg, projlatlon, ref_x, ref_y)
    return [ref_x, ref_y]<jupyter_output><empty_output><jupyter_text>Call geolocation to gather the aoi_coords:<jupyter_code>try:
    aoi_coords = [geolocation(aoi.x1, aoi.y1, geotrans, latlon=False), geolocation(aoi.x2, aoi.y2, geotrans, latlon=False)]
    print(f"aoi_coords in EPSG {utm}: {aoi_coords}")
except TypeError:
    print('TypeError')
    display(Markdown(f'<text style=color:red>This error occurs if an AOI was not selected.</text>'))
    display(Markdown(f'<text style=color:red>Note that the square tool icon in the AOI selector menu is <b>NOT</b> the selection tool. It is the zoom tool.</text>'))
    display(Markdown(f'<text style=color:red>Read the tips above the AOI selector carefully.</text>'))<jupyter_output>aoi_coords in EPSG 32648: [[607882.2025447697, 1724121.6845690133], [768079.9168304839, 1591726.8793742082]]<jupyter_text>Collect the paths to the tiffs:<jupyter_code>tiff_paths = get_tiff_paths(paths)
#print_tiff_paths(tiff_paths)<jupyter_output><empty_output><jupyter_text>Create a subdirectory in which to store the subset tiffs:<jupyter_code># subset_tiffs
print("Choose a directory name in which to store the subset geotiffs.")
print("Note: this will sit alongside the directory containing your pre-subset geotiffs.")
while True:
    sub_name = input()
    if sub_name == "":
        print("Please enter a valid directory name")
        continue
    else:
        break<jupyter_output>Choose a directory name in which to store the subset geotiffs.
Note: this will sit alongside the directory containing your pre-subset geotiffs.
subset_tiffs<jupyter_text>Subset the tiffs and move them from the individual product directories into their own directory, /tiffs:<jupyter_code>tiff_paths = get_tiff_paths(paths)
for p in tiff_paths:
    print(p)
subset_dir = f"{analysis_dir}/{sub_name}/"
new_directory(subset_dir)
for i, tiff_path in enumerate(tiff_paths):
    print(tiff_path)
    date = os.path.basename(date_from_product_name(tiff_path)).split('T')[0]
    polar = get_polarity_from_path(tiff_path)
    print(f"\nProduct #{i+1}:")
    gdal_command = f"gdal_translate -projwin {aoi_coords[0][0]} {aoi_coords[0][1]} {aoi_coords[1][0]} {aoi_coords[1][1]} -projwin_srs 'EPSG:{utm}' -co \"COMPRESS=DEFLATE\" -a_nodata 0 {tiff_path} {subset_dir}{date}_{polar}.tiff"
    print(f"Calling the command: {gdal_command}")
    !{gdal_command}<jupyter_output>Created: /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/subset_tiffs/
/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180619T110424_DVP_RTC30_G_gpufem_3794_VH.tif

Product #1:
Calling the command: gdal_translate -projwin 607882.2025447697 1724121.6845690133 768079.9168304839 1591726.8793742082 -projwin_srs 'EPSG:32648' -co "COMPRESS=DEFLATE" -a_nodata 0 /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180619T110424_DVP_RTC30_G_gpufem_3794_VH.tif /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/subset_tiffs/20180619_VH.tiff
Input file size is 9604, 7299
0...10...20...30...40...50...60...70...80...90...100 - done.
/home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/tiffs/newS1A_IW_20180619T110424_DVP_RTC30_G_gpufem_3794_VV.tif

Product #2:
Calling the command: gdal_translate -projwin 607882.2025447697 1724121.6845690133 768079.9168304839 1591726.8793742082 -projwin_srs 'EPSG:326[...]<jupyter_text>Delete any subset tifs that are filled with NaNs and contain no data.<jupyter_code>subset_paths = f"{subset_dir}*.tif*"
tiff_paths = get_tiff_paths(subset_paths)
remove_nan_filled_tifs(subset_dir, tiff_paths)<jupyter_output>GeoTiffs Examined: 20
GeoTiffs Removed:  0<jupyter_text>Sometimes, when using gdal translate to subset a stack of images, there will be slight differences in sizes of the resulting images, off by a single pixel in either direction. The following code checks the newly subset stack for this problem, and if found, it re-subsets all the images to the size of the smallest image in the stack.<jupyter_code>corners = [gdal.Info(tiff_paths[0], format='json')['cornerCoordinates']['upperLeft'],
                 gdal.Info(tiff_paths[0], format='json')['cornerCoordinates']['lowerRight']]
sizes = list()
for p in tiff_paths:
    info = gdal.Info(p, format='json')
    sizes.append((info['size'][0], info['size'][1]))
    corners = [[max(corners[0][0], gdal.Info(tiff_paths[0], format='json')['cornerCoordinates']['upperLeft'][0]),
               min(corners[0][1], gdal.Info(tiff_paths[0], format='json')['cornerCoordinates']['upperLeft'][1])],
               [min(corners[1][0], gdal.Info(tiff_paths[0], format='json')['cornerCoordinates']['lowerRight'][0]),
               max(corners[1][1], gdal.Info(tiff_paths[0], format='json')['cornerCoordinates']['lowerRight'][1])]]
if len(sizes) != len(set(sizes)):
    print(set(sizes))
    print(corners)
    temp_path = f"{analysis_dir}/temp_subsetting"
    if not os.path.exists(temp_path):
        os.rename(subset_dir, temp_path)
        os.mkdir(subset_dir)
        tiff_paths = glob.glob(f"{temp_path}/*.tif*")

    for i, tiff_path in enumerate(tiff_paths):
        print(f"\nProduct #{i+1}:")
        gdal_command = f"gdal_translate -projwin {corners[0][0]} {corners[0][1]} {corners[1][0]} {corners[1][1]} -projwin_srs 'EPSG:{utm}' -co \"COMPRESS=DEFLATE\" -a_nodata 0 {tiff_path} {subset_dir}{os.path.basename(tiff_path)}"
        print(f"Calling the command: {gdal_command}")
        !{gdal_command}
    shutil.rmtree(temp_path)<jupyter_output>{(5340, 4413)}
[[607860.0, 1724130.0], [768060.0, 1591740.0]]

Product #1:
Calling the command: gdal_translate -projwin 607860.0 1724130.0 768060.0 1591740.0 -projwin_srs 'EPSG:32648' -co "COMPRESS=DEFLATE" -a_nodata 0 /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/temp_subsetting/20180717_VH.tiff /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/subset_tiffs/20180717_VH.tiff
Input file size is 5340, 4413
0...10...20...30...40...50...60...70...80...90...100 - done.

Product #2:
Calling the command: gdal_translate -projwin 607860.0 1724130.0 768060.0 1591740.0 -projwin_srs 'EPSG:32648' -co "COMPRESS=DEFLATE" -a_nodata 0 /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/temp_subsetting/20180806_VV.tiff /home/jovyan/notebooks/ASF/GEOS_657_Project/gamma0_100km_buffer/subset_tiffs/20180806_VV.tiff
Input file size is 5340, 4413
0...10...20...30...40...50...60...70...80...90...100 - done.

Product #3:
Calling the command: gdal_translate -projwin 60[...]<jupyter_text>Decide whether or not to cleanup the original tiffs:<jupyter_code>cleanup = select_parameter(["Save original tiffs", "Delete original tiffs"], '')
cleanup
if cleanup.value == 'Delete original tiffs':
    shutil.rmtree(tiff_dir)<jupyter_output><empty_output><jupyter_text>Print the path to your subset directory:<jupyter_code>print(subset_dir[:-1])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11600055702548391, 'fraction_numerical': 0.07025483915889152, 'mean_word_length': 5.6649651972157775, 'pattern_counts': {'":': 1, '<': 87, '<?xml version=': 0, '>': 86, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '411176', 'license': 'MIT', 'n_tokens_mistral': 6076, 'n_tokens_neox': 5405, 'n_words': 1055, 'path': '2.Subset_Data_Stack.ipynb', 'repo_name': 'biplovbhandari/geos_657_project', 'size': '14372'}