Python script to process Sentinel-2 data
Jump to navigation
Jump to search
Python script to process Sentinel-2 data
Part 1
from qgis.core import (
QgsRasterLayer,
QgsProcessingFeedback,
QgsRectangle
)
from qgis import processing
from osgeo import gdal
import subprocess
import os
import math
os.environ["GDAL_PAM_ENABLED"] = "NO"
# -----------------------------
# USER SETTINGS
# -----------------------------
path = '/mnt/Windows/'
country = 'Alabama'
part = '16R_20250101-20251231'
################################################# IMPORTANT #############################################
# Cuts Sentinel-2 chunks into FlightGear tiles
# Eliminates water only tiles
# Reclassifies Sentenial-2 landcover types to match FlightGear
# Cleans up any urban area road artifacts
# Beginning filename of the NLCD from the source needs to be
# 'nxxxxxxxxxxxxxxxx.tif'
# Also note that "/data/Sentinel-2/" is hardcoded into the path as well
# I needed the existing path structure to process it by state or country.
# Example beginning Sentinel-2 source:
# 16R_20250101-20251231.tif
# Full input path and filename to source example is:
# G:/Scenery/ws3.0/Alabama/data/Sentinel-2/16R_20250101-20251231.tif
# Final output path to processed files is:
# G:/Scenery/ws3.0/Alabama/data/Sentinel-2/tiles_16R_20250101-20251231/N30W086.tif * tile chunk
##########################################################################################################
################################################################### Step 1: Warp land cover 4326 ############################################################
input_raster = path + country + '/data/Sentinel-2/' + part + '.tif'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Land_Cover_4326.tiff'
processing.run(
"gdal:warpreproject",
{
'INPUT':input_raster,
'SOURCE_CRS': QgsCoordinateReferenceSystem('EPSG:5070'),
'TARGET_CRS': QgsCoordinateReferenceSystem('EPSG:4326'),
'RESAMPLING':0,
'NODATA':None,
'TARGET_RESOLUTION':None,
'OPTIONS':'',
'DATA_TYPE':1, # GDAL GDT_Byte
'TARGET_EXTENT':None,
'TARGET_EXTENT_CRS':None,
'MULTITHREADING':False,
'EXTRA':'',
'OUTPUT':output_path
}
)
print("Step 1: Warp land_cover completed. (Sentinel-2_" + country + part + "_Land_Cover_4326)")
############################################################## Step 2: Replace urban and clutter with grass ###################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Land_Cover_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Grass-Only_4326.tiff'
expression = (
'(A == 0) * 44 + '
'(A == 1) * 41 + '
'(A == 2) * 23 + '
'(A == 3) * 44 + '
'(A == 4) * 35 + '
'(A == 5) * 21 + '
'(A == 6) * 44 + '
'(A == 7) * 26 + '
'(A == 8) * 31 + '
'(A == 9) * 34 + '
'(A == 10) * 26 + '
'(A == 11) * 26'
)
command = [
'gdal_calc.py',
'-A', input_path,
'--outfile', output_path,
'--calc', expression,
'--type', 'Byte'
]
subprocess.run(command)
#Windows requires
#subprocess.run(command, shell=True)
print("Step 2: Replace urban and clutter with grass completed. (Sentinel-2_" + country + part + "_Grass-Only_4326)")
################################################################## Step 3: Reclass urban 2 ###########################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Land_Cover_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Reclassed-Urban_4326.tiff'
expression = (
'(A == 7)*2'
)
command = [
'gdal_calc.py',
'-A', input_path,
'--outfile', output_path,
'--calc', expression,
'--type', 'Byte'
]
subprocess.run(command)
#Windows requires
#subprocess.run(command, shell=True)
print("Step 3: Reclass urban completed. (Sentinel-2_" + country + part + "_Reclassed-Urban_4326)")
################################################################ Step 4: Remove clutter and roads from urban ###########################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Reclassed-Urban_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Urban-Only_4326.tiff'
processing.run(
"grass7:r.neighbors",
{
'input':input_path,
'selection':None,
'method':1, #1=median, 2=mode
'size':7,
'gauss':None,
'quantile':'',
'-c':False,
'-a':False,
'weight':'',
'output':output_path,
'GRASS_REGION_PARAMETER':None,
'GRASS_REGION_CELLSIZE_PARAMETER':0,
'GRASS_RASTER_FORMAT_OPT': 'TYPE=BYTE',
'GRASS_RASTER_FORMAT_META':''
}
)
print("Step 4: Remove clutter and roads from urban. (Sentinel-2_" + country + part + "_Urban-Only_4326)")
############################################################## Step 5: Combine grass only and clean urban #########################################################
input_path_a = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Urban-Only_4326.tiff'
input_path_b = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Grass-Only_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Combined-Clean_4326.tiff'
expression = (
'((A > 0) & (B > 0)) * A + '
'((A <= 0) | (B <= 0)) * B'
)
command = [
'gdal_calc.py',
'-A', input_path_a,
'-B', input_path_b,
'--outfile', output_path,
'--calc', expression,
'--NoDataValue', '0',
'--type', 'Byte'
]
subprocess.run(command)
#Windows requires
#subprocess.run(command, shell=True)
print("Step 5: Combined and clean completed. (Sentinel-2_" + country + part + "_Combined-Clean_4326)")
############################################################## Step 6: Cut into tiles #########################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Combined-Clean_4326.tiff'
output_folder = path + country + '/data/Sentinel-2/tiles_' + part + '/'
tile_size_deg = 1.0
os.makedirs(output_folder, exist_ok=True)
rlayer = QgsRasterLayer(input_path, "input")
if not rlayer.isValid():
raise Exception("Raster failed to load: " + input_path)
ext = rlayer.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
xmin_int = math.floor(xmin)
xmax_int = math.ceil(xmax)
ymin_int = math.floor(ymin)
ymax_int = math.ceil(ymax)
feedback = QgsProcessingFeedback()
for lat in range(ymin_int, ymax_int - 1):
for lon in range(xmin_int, xmax_int - 1):
tile_xmin = lon
tile_xmax = lon + tile_size_deg
tile_ymin = lat
tile_ymax = lat + tile_size_deg
if tile_xmax <= xmin or tile_xmin >= xmax:
continue
if tile_ymax <= ymin or tile_ymin >= ymax:
continue
tile_extent = f"{tile_xmin},{tile_xmax},{tile_ymin},{tile_ymax}"
ns = "N" if (lat + 1) >= 0 else "S"
ew = "E" if lon >= 0 else "W"
out_name = f"{ns}{abs(lat+1):02d}{ew}{abs(lon):03d}.tif"
out_path = os.path.join(output_folder, out_name)
tmp_tile = out_path.replace(".tif", "_tmp.tif")
processing.run(
"gdal:cliprasterbyextent",
{
'INPUT': input_path,
'PROJWIN': tile_extent,
'NODATA': 0,
'OPTIONS': '',
'DATA_TYPE': 0,
'OUTPUT': tmp_tile
},
feedback=feedback
)
if not os.path.exists(tmp_tile) or os.path.getsize(tmp_tile) < 50:
if os.path.exists(tmp_tile):
os.remove(tmp_tile)
print("Skipped (empty tile):", out_name)
continue
ds = gdal.Open(tmp_tile)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
# If every pixel is zero → water only
if (arr == 0).all():
os.remove(tmp_tile)
print("Skipped (all water):", out_name)
continue
os.rename(tmp_tile, out_path)
print("Created:", out_path)
print("Tiling complete.")
Part 2
from qgis.core import QgsApplication, QgsCoordinateTransform, QgsProject, QgsRasterLayer, QgsCoordinateReferenceSystem, QgsProcessingException, QgsRasterBlock, QgsRectangle
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis import processing
from processing.core.Processing import Processing
from osgeo import gdal, osr, ogr
import os
import numpy
import numpy as np
import subprocess
# Define input layer names, change according to your file names
path = '/mnt/Windows/';
country = 'Alabama';
part = 'N30W086';
#ratio of upscale to smoothing
percentage = 3
size = 11
################################## IMPORTANT ###############################
# Beginning filename of the NLCD from the source needs to be
# 'nxxxxxxxxxxxxxxxx.tif'
# Also note that "/data/Sentinel-2/" is hardcoded into the path as well
# I needed the existing path structure to process it by state or country.
# Example beginning Sentinel-2 source:
# N30W086.tif
# Full input path and filename to source example is:
# G:/Scenery/ws3.0/Alabama/data/Sentinel-2/N30W086.tif
# Final output path to last processed file is:
# G:/Scenery/ws3.0/Alabama/data/Sentinel-2/Sentinel-2_AlabamaN30W086_Smoothed-HD-Compressed_4326.tiff
############################################################################
################################################################## Step 1: Upsample to HD #################################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Combined-Clean_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Combined-Clean-HD_4326.tiff'
# Open the original raster
ds = gdal.Open(input_path)
gt = ds.GetGeoTransform()
# Extract the original resolution
original_xRes = gt[1]
original_yRes = abs(gt[5])
# Define the percentage to resize by (e.g., 0.50 for 50%, 2.0 for 200%)
#defined at the begining of the script
#percentage = 8.0 # Adjust this as needed
# Calculate the new resolution
new_xRes = original_xRes / percentage
new_yRes = original_yRes / percentage
# Perform the warp (resampling)
gdal.Warp(
output_path,
input_path,
xRes=new_xRes,
yRes=new_yRes,
outputType=gdal.GDT_Byte # Set the output type to uint8
)
print("Step 1: Upsample to. (Sentinel-2_" + country + part + "_Combined-Clean-HD_4326)")
################################################################## Step 2: Smooth all features in original dataset #################################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Combined-Clean-HD_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Smoothed-HD_4326.tiff'
processing.run(
"grass7:r.neighbors",
{
'input':input_path,
'selection':None,
'method':1, #1=median, 3=mode
'size':size,
'gauss':None,
'quantile':'',
'-c':False,
'-a':False,
'weight':'',
'output':output_path,
'GRASS_REGION_CELLSIZE_PARAMETER':0,
'GRASS_RASTER_FORMAT_OPT': 'TYPE=BYTE',
'GRASS_RASTER_FORMAT_META':''
}
)
print("Step 2: Smooth all features in original dataset completed. (Sentinel-2_" + country + part + "_Smoothed-HD_4326)")
################################################################## Step 3: Convert to 8Bit #################################################################
input_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Smoothed-HD_4326.tiff'
output_path = path + country + '/data/Sentinel-2/Sentinel-2_' + country + part + '_Smoothed-HD-Compressed_4326.tiff'
processing.run(
"gdal:translate",
{
'INPUT':input_path,
'TARGET_CRS':None,
'NODATA':0,
'COPY_SUBDATASETS':False,
'OPTIONS':'',
'EXTRA':'',
'DATA_TYPE':1,
'OUTPUT':output_path,
'OPTIONS': 'COMPRESS=LZW'
}
)
result_bit_conversion_layer = QgsRasterLayer(output_path, 'Sentinel-2_' + country + part + '_Smoothed-HD-Compressed_4326')
QgsProject.instance().addMapLayer(result_bit_conversion_layer)
print("Step 3: Convert to Compressed. (Sentinel-2_" + country + part + "_Smoothed-HD-Compressed_4326)")