<jupyter_start><jupyter_text>Landsat Imagery Preprocessing: Clip raster and fill gaps<jupyter_code>import numpy as np
import math
import scipy
import pandas as pd
import PIL
import gdal
import matplotlib.pyplot as plt
plt.style.use('ggplot')
import sys, os
from pathlib import Path
import time
import random
import collections, functools, operator
import csv
import subprocess
import datetime

from osgeo import gdal,osr
from gdalconst import *
import subprocess
from osgeo.gdalconst import GA_Update<jupyter_output><empty_output><jupyter_text>Clip raster function<jupyter_code>#=============================[ CLIP RASTER ]=================================#

#Clips a selected input raster to have the exact same projection, extent, and resolution as another
# USAGE: python clip_to_raster.py src clipsrc dst
# INPUTS: src - the raster to reproject
#         clipsrc - the raster to match
# OUTPUT: dst - the output file


# Important: If the pixels in src and clipsrc do not line up exactly, the pixels in SRS will be resampled
# To preserve the extents, even if they have the exact same resolution

def clipRaster(src, clipsrc, dst):
    # Parameters
    resampleOptions = ('average','near','bilinear','cubic','cubicspline','lanczos','mode','max','min','med','Q1','Q3')
    resample = 'cubic'
    tr = ''
    ts = ''
    dstnodata = ''
    overwrite = False


    # Test if there is a problem opening the input data
    ds1 = gdal.Open(src, GA_ReadOnly)
    if ds1 is None:
        print('Could not open ' + src)
        sys.exit(1)
        
    ds2 = gdal.Open(clipsrc, GA_ReadOnly)
    if ds2 is None:
        print ('Could not open ' + clipsrc)
        sys.exit(1)
        
    # Read the georeferecing info on the raster_to_match
    transform1 = ds1.GetGeoTransform()
    transform2 = ds2.GetGeoTransform()
    wkt = ds2.GetProjection()
    rows = ds2.RasterYSize
    cols = ds2.RasterXSize

    pixelWidth = transform1[1]
    pixelHeight = transform1[5]

    minX = transform2[0]
    maxY = transform2[3]
    maxX = minX + (cols * pixelWidth)
    minY = maxY + (rows * pixelHeight)

    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    proj4 = srs.ExportToProj4()
        
    # Build Warp Options for gdalwarp based on the program inputs
    warpOptions = gdal.WarpOptions(format='GTiff', xRes=pixelWidth, yRes=pixelHeight,
                    srcNodata=0, dstNodata='nan', width=0, height=0, outputBounds=[minX, minY, maxX, maxY],
                    resampleAlg=resample, multithread=True, dstSRS=proj4, cropToCutline=True)

    # Reproject the input_raster to match the raster_to_match
    outBand = gdal.Warp(dst, src, options=warpOptions)
    outBand = None<jupyter_output><empty_output><jupyter_text>Examples with data from USGS<jupyter_code>listdir_src = os.getcwd() + '/drive/My Drive/TFG/Ciudad Real/test/EO-1 ALI/'
listdir_src = [os.path.join(listdir_src, f) for f in sorted(os.listdir(listdir_src))]
clipsrc = os.getcwd() + '/drive/My Drive/TFG/Ciudad Real/test/EO-1 Hyperion/'
clipsrc = os.getcwd() + '/drive/My Drive/TFG/Ciudad Real/test/EO-1 Hyperion/' + os.listdir(clipsrc)[0]
for src in listdir_src:
    clipRaster(src, clipsrc, src)


# Directories
listdir_src = os.getcwd() + '/drive/My Drive/TFG/Karnataka_20scenes/Landsat 7/'
listdir_dst = os.getcwd() + '/drive/My Drive/TFG/Karnataka_20scenes/Landsat 7 Processed/'

listdir_dst = [os.path.join(listdir_dst, f) for f in sorted(os.listdir(listdir_src))]
listdir_src = [os.path.join(listdir_src, f) for f in sorted(os.listdir(listdir_src))]

listfile_clipsrc = os.getcwd() + '/drive/My Drive/TFG/Karnataka_20scenes/EO-1 Hyperion/'
listfile_clipsrc = [os.path.join(listfile_clipsrc, f) for f in sorted(os.listdir(listfile_clipsrc))]
listfile_clipsrc = [f + '/' + f.split('/')[-1][:-5] + '/' + f.split('/')[-1][:-4]+'B008_L1GST.TIF' for f in listfile_clipsrc]


driver = gdal.GetDriverByName('GTiff')



#=============[ PROCESSING OF EACH BAND: fill gaps & cropping ]================#
#FOR USGS
for dir_src, file_clipsrc, dir_dst in zip(listdir_src, listfile_clipsrc, listdir_dst):
    Path(dir_dst).mkdir(parents=True, exist_ok=True)
    for bandName in os.listdir(dir_src):
        if 'band' in bandName: # sustituir esto por alguna funcion que busque en el metadata si es banda o no. Esto solo funciona para landsat
            file_src = os.path.join(dir_src, bandName)
            file_dst = os.path.join(dir_dst, bandName)

            src_ras = gdal.Open(file_src, GA_Update)
            ras = driver.CreateCopy(file_dst, src_ras, strict=0)
            src_ras = None
            
            band = ras.GetRasterBand(1)
            gdal.FillNodata(targetBand=band, maskBand=None, maxSearchDist = 20, smoothingIterations = 0) # Fills gaps
            ras = None # unlink the file object and save the results. For gdal to write something to disk you need to flush/close the dataset.
            clipRaster(file_dst, file_clipsrc, file_dst) # Clips raster

# For increasing performance I could first clip the raster with a margin for performing fillNoData, then do fillNoData and then clip again to match Hyperion
# Directories
listdir_src = os.getcwd() + '/drive/My Drive/TFG/Ciudad Real/EO-1 ALI/'
listdir_dst = os.getcwd() + '/drive/My Drive/TFG/Ciudad Real/EO-1 ALI Processed/'

listdir_dst = [os.path.join(listdir_dst, f, f[:-5]) for f in sorted(os.listdir(listdir_src))]
listdir_src = [os.path.join(listdir_src, f, f[:-5]) for f in sorted(os.listdir(listdir_src))]

listfile_clipsrc = os.getcwd() + '/drive/My Drive/TFG/Ciudad Real/EO-1 Hyperion/'
listfile_clipsrc = [os.path.join(listfile_clipsrc, f, f[:-5], f[:-4]+'B008_L1GST.TIF') for f in sorted(os.listdir(listfile_clipsrc))]

driver = gdal.GetDriverByName('GTiff')

bandsToGet = ['B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B10']


#=============[ PROCESSING OF EACH BAND: cropping ]================#
#FOR USGS
for dir_src, file_clipsrc, dir_dst in zip(listdir_src, listfile_clipsrc, listdir_dst):
    Path(dir_dst).mkdir(parents=True, exist_ok=True)
    for bandName in os.listdir(dir_src):
        if any(b in bandName for b in bandsToGet):
            file_src = os.path.join(dir_src, bandName)
            file_dst = os.path.join(dir_dst, bandName)
            clipRaster(file_src, file_clipsrc, file_dst) # Clips raster<jupyter_output><empty_output><jupyter_text>Example with data from Google Earth Engine (unfinished)<jupyter_code>#FOR EARTH ENGINE

dir_dst += 'LE07_144052_20140627.tif'
clipRaster(dir_src, dir_clipsrc, dir_dst) # Clips raster
#tropt = gdal.TranslateOptions(noData=0., bandList=[1,2,3,4,5,6], format="GTiff")
#gdal.Translate(dir_dst, dir_dst, options=tropt)
directory = 'drive/MyDrive/TFG/Karnataka_20scenes/Landsat 7 Processed/Not masked'
#%cd drive/MyDrive/TFG/Karnataka_20scenes/Landsat 7 Processed/Masked
#!gdal_translate -of GTiff -a_nodata 0 LE07_144052_20140627.tif LE07_144052_20140627_wNoData.tif

dir_dst = os.getcwd() + '/drive/My Drive/TFG/Karnataka_20scenes/Landsat 7 Processed/Not masked/LE07_144052_20140627.tif'


ras = gdal.Open(dir_dst, GA_Update)
for i in range(1,7):
    band = ras.GetRasterBand(i)
    bandarr = band.ReadAsArray()
    gdal.FillNodata(targetBand=band, maskBand=None, maxSearchDist = 40, smoothingIterations = 0) # Fills gaps
    band=None
ras = None



im = gdal.Open(dir_dst, GA_Update)
r = im.GetRasterBand(3).ReadAsArray()
g = im.GetRasterBand(2).ReadAsArray()
b = im.GetRasterBand(1).ReadAsArray()
plt.hist(r.flatten(), bins=200)
plt.show()
rgb = normalize_bands([r,g,b])
plt.hist(rgb[0].flatten(), bins=200)
plt.show()
viz_raster(rgb[0], rgb[1], rgb[2])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13000513347022588, 'fraction_numerical': 0.023485626283367558, 'mean_word_length': 3.7547284929835265, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '646534', 'license': 'MIT', 'n_tokens_mistral': 3031, 'n_tokens_neox': 2849, 'n_words': 656, 'path': '2.Datasets/Landsat_Imagery_Preprocessing.ipynb', 'repo_name': 'Rojas-D/thesis', 'size': '7792'}