ricampiona i file jp2 originali di Sentinel 2 in geotiff a 10 m
import os
import rasterio
import numpy as np
from rasterio.windows import Window
from pathlib import Path
def create_tiles_from_geotiffs(input_files, output_dir, tile_size_km=2.5):
"""
Cut multiple geotiffs into tiles and stack them as bands in new geotiffs.
Args:
input_files: List of paths to input geotiffs (must cover same area)
output_dir: Directory to save output tiles
tile_size_km: Size of tiles in kilometers
"""
# Create output directory if it doesn't exist
os.makedirs(output_dir, exist_ok=True)
# Open first file to get metadata
with rasterio.open(input_files[0]) as src:
# Get spatial reference metadata
crs = src.crs
transform = src.transform
# Calculate tile size in pixels
# Get resolution in meters per pixel
res_x = abs(transform.a) # x resolution in map units/pixel
res_y = abs(transform.e) # y resolution in map units/pixel
# Calculate tile size in pixels
tile_size_x = int(tile_size_km * 1000 / res_x)
tile_size_y = int(tile_size_km * 1000 / res_y)
# Get image dimensions
width = src.width
height = src.height
# Calculate number of tiles
num_tiles_x = (width + tile_size_x - 1) // tile_size_x
num_tiles_y = (height + tile_size_y - 1) // tile_size_y
print(f"Input size: {width}x{height} pixels")
print(f"Tile size: {tile_size_x}x{tile_size_y} pixels ({tile_size_km}km x {tile_size_km}km)")
print(f"Creating {num_tiles_x * num_tiles_y} tiles ({num_tiles_x}x{num_tiles_y})")
# Process each tile
for tile_y in range(num_tiles_y):
for tile_x in range(num_tiles_x):
# Calculate window coordinates
x_offset = tile_x * tile_size_x
y_offset = tile_y * tile_size_y
# Adjust for edge tiles
actual_tile_width = min(tile_size_x, width - x_offset)
actual_tile_height = min(tile_size_y, height - y_offset)
# Skip if tile is empty or too small
if actual_tile_width <= 0 or actual_tile_height <= 0:
continue
# Define the window to read
window = Window(x_offset, y_offset, actual_tile_width, actual_tile_height)
# Calculate new geotransform for the tile
new_transform = rasterio.transform.from_origin(
transform.c + transform.a * x_offset, # upper left x
transform.f + transform.e * y_offset, # upper left y
transform.a, # pixel width
transform.e # pixel height
)
# Create empty array for the 11-band tile
tile_data = np.zeros((len(input_files), actual_tile_height, actual_tile_width),
dtype=rasterio.float32)
# Read data from each input file
for band_idx, file_path in enumerate(input_files):
with rasterio.open(file_path) as src:
tile_data[band_idx] = src.read(1, window=window)
# Define output file path
out_path = os.path.join(output_dir, f"tile_x{tile_x}_y{tile_y}.tif")
# Write the tile to a new file
with rasterio.open(
out_path,
'w',
driver='GTiff',
height=actual_tile_height,
width=actual_tile_width,
count=len(input_files),
dtype=tile_data.dtype,
crs=crs,
transform=new_transform,
) as dst:
dst.write(tile_data)
print(f"Created tile: {out_path}")
# Optional: Add metadata about original tile position
with rasterio.open(out_path, 'r+') as dst:
dst.update_tags(
tile_x=tile_x,
tile_y=tile_y,
original_extent=f"{width}x{height}",
tile_size_km=tile_size_km
)
if __name__ == "__main__":
# Directory containing input geotiffs
input_dir = "output_10m"
# List of input geotiff files (11 bands)
input_files = [
os.path.join(input_dir, f"B{i+1}.tif") for i in range(11)
]
# Check if files exist
for f in input_files:
if not os.path.exists(f):
print(f"Warning: File {f} does not exist!")
# Directory to save output tiles
output_dir = "output_tiles"
# Create tiles
create_tiles_from_geotiffs(input_files, output_dir, tile_size_km=2.56)
crea una tile di 2560x2560m 11 bande
import os
import rasterio
from rasterio.windows import Window
import numpy as np
# Parameters
input_folder = "./output_10m" # folder with B1.tif to B11.tif
output_folder = "./output_tiles"
tile_size_m = 2560 # meters
pixel_size = 10 # meters per pixel
tile_size_px = tile_size_m // pixel_size
# Get list of sorted input files (assumes filenames like B1.tif, B2.tif, ..., B11.tif)
input_files = sorted([f for f in os.listdir(input_folder) if f.endswith(".tif")])
input_paths = [os.path.join(input_folder, f) for f in input_files]
# Open all rasters
rasters = [rasterio.open(fp) for fp in input_paths]
# Check that all rasters have same dimensions
width, height = rasters[0].width, rasters[0].height
transform = rasters[0].transform
crs = rasters[0].crs
# Create output folder
os.makedirs(output_folder, exist_ok=True)
tile_id = 0
for y in range(0, height, tile_size_px):
for x in range(0, width, tile_size_px):
if x + tile_size_px > width or y + tile_size_px > height:
continue # Skip partial tiles
# Read corresponding tile from each band
tile_stack = []
for r in rasters:
tile = r.read(1, window=Window(x, y, tile_size_px, tile_size_px))
tile_stack.append(tile)
tile_stack = np.stack(tile_stack) # shape: (bands, height, width)
# Define transform for this tile
tile_transform = rasterio.windows.transform(Window(x, y, tile_size_px, tile_size_px), transform)
# Write output GeoTIFF with multiple bands
out_path = os.path.join(output_folder, f"tile_{tile_id:04d}.tif")
with rasterio.open(
out_path,
"w",
driver="GTiff",
height=tile_size_px,
width=tile_size_px,
count=len(rasters),
dtype=tile_stack.dtype,
crs=crs,
transform=tile_transform,
) as dst:
for i in range(len(rasters)):
dst.write(tile_stack[i], i + 1)
tile_id += 1
# Close all files
for r in rasters:
r.close()
print(f"Done. Created {tile_id} multi-band tiles.")
