<jupyter_start><jupyter_text>Split Sentinel-2 data to tiles<jupyter_code># Imports
import rasterio
from rasterio.plot import show
from rasterio.windows import Window, bounds, from_bounds 
from rasterio.transform import rowcol, xy

from sentinelhub import BBox, CRS
from pyproj import Transformer
from shapely.geometry import Polygon

import h5py
import numpy as np

import time
import json
import glob
import os
from tqdm.notebook import tqdm<jupyter_output><empty_output><jupyter_text>Export a window of Landsat-2 data<jupyter_code># File and folder paths
dirpath = r"data\sentinel2"
search_criteria = "*\\R10m\\transformed_tci.jp2"
q = os.path.join(dirpath, search_criteria)
print(q)
filenames = glob.glob(q)
filenames
files = []
for filename in filenames:
    file = rasterio.open(filename)
    files.append(file)
    
files
# given a window and some bounds (left, bottom, right, top)
# checks if the window in contained inside the bounds
def window_in_bounds(window, bounds):
    w_left, w_bottom, w_right, w_top = window
    b_left, b_bottom, b_right, b_top = bounds
    
    in_horizontal = w_left > b_left and w_right < b_right
    in_vertical = w_bottom > b_bottom and w_top < b_top
    
    return in_horizontal and in_vertical

# Given the metadata file (tileInfo.json) from a landsat-2 raster, 
# computer the wgs84 aligned inner bounds contained inside it
def inner_bounds_from_metadata(metadata):
    data = json.load(open(metadata))    

    crs_name = data["tileGeometry"]["crs"]["properties"]["name"].split(':')[-3]
    crs_number = data["tileGeometry"]["crs"]["properties"]["name"].split(':')[-1]

    coordinates = data["tileGeometry"]["coordinates"][0]
    transformer = Transformer.from_crs(f'{crs_name}:{crs_number}', "epsg:4326")
    coordinates = [list(transformer.transform(point[0], point[1])) for point in coordinates]
    coordinates = [[point[1], point[0]] for point in coordinates]

    polygon = Polygon(coordinates)
    
    x, y = polygon.exterior.coords.xy
    # sort coords and remove duplicate
    x = sorted(x[:-1])
    y = sorted(y[:-1])
    
    min_x, max_x = x[1], x[2]
    min_y, max_y = y[1], y[2]
    
    return (min_x, min_y, max_x, max_y)
# Given a ls2 file, a wgs84 coordinate, a window size in pixels and a radiance value (optional),
# exports small images centered on the coordinate of window_size containing the 
# coordinate and radiance value in the filename
def export_split(ls2_file, center_of_window, window_size, px=None, py=None, radiance=None):
    filename_splits = ls2_file.name.split('\\')

    input_filename = ls2_file.name     
    metadata_filename = ('\\').join(filename_splits[0:3]) + '\\tileInfo.json'

    file_transform = ls2_file.transform

    # center of window (pixels)
    x, y = rowcol(file_transform, center_of_window[0], center_of_window[1], round)

    # bound limits (pixels)
    left_px, bottom_py = x + window_size/2, y - window_size/2
    right_px, top_py =   x - window_size/2, y + window_size/2
    # bound limits (wgs84)
    left, bottom = xy(file_transform, left_px, bottom_py, offset='center')
    right, top = xy(file_transform, right_px, top_py, offset='center')    
    # bounds (wgs84)
    calc_bounds = (left, bottom, right, top)
    inner_file_bounds = inner_bounds_from_metadata(metadata_filename)

    if window_in_bounds(calc_bounds, inner_file_bounds):
        # Create a Window and calculate the transform from the source dataset
        window = from_bounds(left, bottom, right, top, file_transform)
        window_transform = ls2_file.window_transform(window)

        # Create a new cropped raster to write to
        profile = ls2_file.profile
        profile.update({
            'driver': 'PNG',
            'height': window_size,
            'width': window_size,
            'transform': window_transform})

        output_filename = ('\\').join(filename_splits[0:2]) + f'\\splits\\sokoto\\{center_of_window}_{px}x{py}_{radiance}.png'
        with rasterio.open(output_filename, 'w', **profile) as dst:
            # Read the data from the window and write it to the output image
            dst.write(ls2_file.read(window=window))
        return True # Successfully exported
            
    return False # Not in window<jupyter_output><empty_output><jupyter_text>Open blackmarble file and get infos<jupyter_code># Geo stuff
def tileId_to_geo(tile_id):   
    horizontalId = int(tile_id[1:3])
    verticalId = int(tile_id[4:6])
    
    lat = (horizontalId - 18) * 10
    lon = -(verticalId - 9) * 10
    
    return (lat, lon)

# Origin (px and geo): top left
def px_to_geo_h5(px, py, x_width, y_height, tile_geo_width, tile_geo_height, geo_tile):
    ratio_px = px / x_width
    ratio_py = py / y_height
    
    geo_dx = ratio_px * tile_geo_width
    geo_dy = ratio_py * tile_geo_height
    
    return (geo_tile[0] + geo_dx, geo_tile[1] - geo_dy)

def geo_to_px_h5(geo_target, geo_tile, x_width, y_height, tile_geo_width, tile_geo_height):        
    diff_lat  = geo_target[0] - geo_tile[0] # lat of tile will always be smaller
    diff_long = geo_tile[1] - geo_target[1] # long of tile always larger
    
    if not (diff_lat >= 0 and diff_lat <= tile_geo_height and diff_long >= 0 and diff_long <= tile_geo_width):
        raise Exception(f'Target geo point not in provided tile: {geo_target} not in {geo_tile}, {geo_tile[0] + tile_geo_width}, {geo_tile[1] - tile_geo_height} ({tile_geo_width}°x{tile_geo_height}°)')
    
    ratio_px = diff_lat / tile_geo_height
    ratio_py = diff_long / tile_geo_width
    
    return (int(ratio_px * x_width), int(ratio_py * y_height))
# Blackmarble H5 File
def geo_tl(h5file):
    splits = h5file.split('.')
    
    tile_id = splits[2]
    geo = tileId_to_geo(tile_id)                              
    
    return geo

def export_data_h5(h5file):
    name='DNB_At_Sensor_Radiance_500m'
    path='/HDFEOS/GRIDS/VNP_Grid_DNB/Data Fields/'
    data_dir='data/blackmarble/'
    
    with h5py.File(data_dir + h5file, 'r') as f:
        data = np.array(f[path + name][:,:])
        
        return data<jupyter_output><empty_output><jupyter_text>Traverse h5 file and for each pixel, export the corresponding Landsat-2 window (if it exists)<jupyter_code># Time estimation:
def time_estimation(xlim, ylim):
    import datetime
    xmin, xmax = xlim.start, xlim.stop
    ymin, ymax = ylim.start, ylim.stop
    time_per_pxl_file = 1866 / 59950 # Empirical value (1866s for 59950 pixels)
    pxls = (xmax-xmin) * (ymax-ymin)
    secs = time_per_pxl_file * pxls
    time = datetime.timedelta(seconds=secs)
    
    print(f'Estimated {time}')
h5file = 'VNP46A1.A2021165.h18v07.001.2021166080518.h5'

tile_width_deg = 10
split_width_px = 50

data = export_data_h5(h5file)
exported_bitmap = np.zeros((2400, 2400), dtype=bool)

geo_tile = geo_tl(h5file)

print(f'Nightlight file: Top Left={geo_tile}, Bottom Right=({geo_tile[0]+tile_width_deg}, {geo_tile[1]-tile_width_deg})')
print()

for file in files:
    print(file)
    left, bottom, right, top = file.bounds[0], file.bounds[1], file.bounds[2], file.bounds[3]
    
    print(f'Degrees: Left: {left}, Bottom: {bottom}, Right: {right}, Top: {top}')
    
    left_px, bottom_px = geo_to_px_h5((left, bottom), geo_tile, 2400, 2400, tile_width_deg, tile_width_deg)
    right_px, top_px = geo_to_px_h5((right, top), geo_tile, 2400, 2400, tile_width_deg, tile_width_deg)
    
    print(f'Pixels: Left: {left_px}, Bottom: {bottom_px}, Right: {right_px}, Top: {top_px}')
    print()
    
    xlim = range(left_px, right_px)
    ylim = range(top_px, bottom_px)
    
    count_exp = 0
    count_oob = 0

    for px in tqdm(xlim):
        for py in ylim:
            if (exported_bitmap[px, py] == True):
                continue
            # for some reason, x and y coordinates are flipped in data
            radiance = data[py, px]
            lat, lon = px_to_geo_h5(px, py, data.shape[0], data.shape[1], tile_width_deg, tile_width_deg, geo_tile)        
            exported = export_split(file, center_of_window = (lat, lon), window_size = split_width_px, px=px, py=py, radiance = radiance)
            
            if exported:
                exported_bitmap[px, py] = 1<jupyter_output>Nightlight file: Top Left=(0, 20), Bottom Right=(10, 10)

<open DatasetReader name='data\sentinel2\31PFP,2021-06-14,0\R10m\transformed_tci.jp2' mode='r'>
Degrees: Left: 3.9174820064733487, Bottom: 11.665122014929599, Right: 4.931722160292178, Top: 12.662956415002258
Pixels: Left: 940, Bottom: 2000, Right: 1183, Top: 1760
{'fraction_non_alphanumeric': 0.1130322277111268, 'fraction_numerical': 0.03363914373088685, 'mean_word_length': 3.596216216216216, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, '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': '432493', 'license': 'MIT', 'n_tokens_mistral': 3205, 'n_tokens_neox': 3026, 'n_words': 738, 'path': 'jupyter-notebooks/3 - split_tile.ipynb', 'repo_name': 'matt989253/heig-TerraGazer', 'size': '8502'}