Python script to process Sentinel-2 data
Jump to navigation
Jump to search
Python script to process Sentinel-2 data
This is a slightly different process that that of the NLCD. Sentinel-2 data does not have any extra tree canopy coverage to merge in but this script still cuts out some superfluous urban clutter. Sentinel-2 data uses a two stage workflow. Step 1 cleans up the urban clutter and then splits the larger Sentinel-2 chunks of date to 1x1 tiles to help facilitate Step 2, which up-samples and smoothes the boundaries.
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
This is the pretty much the same final steps that the NLCD uses to up-sample and smooth the boundaries.
NOTE: the current genVPB.py script (osgdem) will not generate the final scenery water/land edge to this high of a resolution. It is set to a default tile size of 256x256, "--tile-image-size 256". This resolution is much higher than that. To achieve the same base resolution in the built scenery, osgdem needs a switch to "--tile-image-size 1024" to realize the higher resolution. However, when you apply the water mask flags --generate-water-raster --shrink-water 4 --coastline [path_to_coastal_shapefiles] using genVPB.py the cut mask will achieve the higher resolution.
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 = 8
size = 11
#ratio of 11/15 upscale to smoothing is the next sweet spot, but starts to move boundary positions and starts to eliminate smaller land class features for not a ton of resolution gain.
################################## 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 beginning 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)")