import os
import re
import json
import shapely
import rasterio
import datetime
import dateutil
import pystac_client
import multiprocessing
from smosaic.smosaic_clip_raster import clip_raster
from smosaic.smosaic_collection_get_data import collection_get_data
from smosaic.smosaic_collection_query import collection_query
from smosaic.smosaic_count_pixels import count_pixels
from smosaic.smosaic_filter_scenes import filter_scenes
from smosaic.smosaic_fix_baseline_number import fix_baseline_number
from smosaic.smosaic_generate_cog import generate_cog
from smosaic.smosaic_get_dataset_extents import get_dataset_extents
from smosaic.smosaic_grid_crop import clip_from_grid
from smosaic.smosaic_merge_scene import merge_scene, merge_scene_provenance_cloud
from smosaic.smosaic_merge_tifs import merge_tifs
from smosaic.smosaic_reproject_tif import reproject_tifs
from smosaic.smosaic_spectral_indices import calculate_spectral_indices
from smosaic.smosaic_utils import add_days_to_date, add_months_to_date, clean_dir, days_between_dates, get_all_cloud_configs, load_jsons, map_band_name
[docs]
def mosaic(name, data_dir, stac_url, collection, output_dir, start_year, start_month, start_day, mosaic_method, stac_source="bdc", grid_crop=False, bands=None, reference_date=None, duration_days=None, end_year=None, end_month=None, end_day=None, duration_months=None, geom=None, grid=None, tile_id=None, bbox=None, profile=None, projection_output=4326):
"""
Create Analysis-Ready Data (ARD) seamless satellite image mosaics directly from STAC catalogs.
This function generates temporal composites by applying specified mosaic composition functions over
defined time periods and spatial extents. The system connects to STAC catalogs from the
Brazil Data Cube (BDC), Digital Earth Africa, Swiss Data Cube, and Microsoft Planetary Computer.
Args:
name (str): Descriptive identifier for the output mosaic (used in file naming).
data_dir (str): Directory for storing intermediate data and processing artifacts.
stac_url (str): STAC API endpoint URL (e.g., "https://data.inpe.br/bdc/stac/v1").
collection (str): STAC collection identifier (e.g., "S2_L2A-1").
output_dir (str): Destination directory for final mosaic products.
start_year (int): Initial year for temporal filtering (YYYY format).
start_month (int): Initial month for temporal filtering (1-12).
start_day (int): Initial day for temporal filtering (1-31).
mosaic_method (str): Mosaic composition function. Supported methods include:
- "lcf": Least Cloud-cover First - order by the least cloud-cover.
- "chrono": Chronological - order chronologically.
- "ctd": Closest to Date - order by the closest image to reference date.
stac_source (str, optional): The source STAC catalog ecosystem. Supported sources include:
- "bdc": Brazil Data Cube (default)
- "planetary-computer": Microsoft Planetary Computer
- "digitalearth-africa": Digital Earth Africa
- "swissdatacube": Swiss Data Cube
Defaults to "bdc".
grid_crop (bool, optional): Enable cropping to grid tile boundaries. Defaults to False.
bands (list, optional): Spectral bands to include (e.g., ["B02","B03","B04","B8A"]).
reference_date (str, optional): Reference date for the Closest to Date composition function.
(format: 'YYYY-MM-DD').
duration_days (int, optional): Temporal window length in days from start date.
end_year (int, optional): Explicit end year for temporal range (YYYY format). Defaults to None.
end_month (int, optional): Explicit end month for temporal range (1-12). Defaults to None.
end_day (int, optional): Explicit end day for temporal range (1-31). Defaults to None.
duration_months (int, optional): Temporal window length in months from start date.
geom (str/dict, optional): GeoJSON geometry defining Area of Interest (AOI).
Defaults to None.
grid (str, optional): Grid system identifier. Supported grids include:
- "BDC_SM_V2": Brazil Data Cube Small Grid.
Defaults to None.
tile_id (str, optional): Grid tile identifier (e.g., 6-digit code "020019").
Requires grid parameter. Defaults to None.
bbox (list/tuple, optional): Bounding box coordinates [min_lon, min_lat, max_lon, max_lat].
profile (dict, optional): Profile band and spectral indices selector.
Overrides bands parameters (e.g., "urban_analysis" or "crop_condition").
projection_output (int/str, optional): Output coordinate reference system. Options:
- EPSG codes: 4326 (WGS84), 5880 (SIRGAS 2000 Brazil Polyconic)
- BDC codes: "BDC" (Brazil Data Cube Standard Grid projection)
Defaults to 4326.
Example:
>>> import os
>>> from smosaic import mosaic
>>>
>>> stac_url = "https://data.inpe.br/bdc/stac/v1"
>>>
>>> result = mosaic(
... name="MT",
... data_dir=os.path.abspath(""),
... stac_url=stac_url,
... collection="S2_L2A-1",
... grid="BDC_SM_V2",
... tile_id="020019",
... projection_output="BDC",
... output_dir=os.path.join("output"),
... mosaic_method="lcf",
... start_year=2026,
... start_month=1,
... start_day=1,
... duration_days=16,
... bands=["B02","B03","B04"]
... )
Returns:
str: Path to the generated mosaic file.
"""
clean_dir(data_dir)
stac = pystac_client.Client.open(stac_url)
if (stac_source == "bdc" and collection == "S2_L1C_BUNDLE-1" or
stac_source == "bdc" and collection == "S2_L2A-1" or
stac_source == "planetary-computer" and collection == "sentinel-2-l2a" or
stac_source == "digitalearth-africa" and collection == "s2_l2a" or
stac_source == "swissdatacube" and collection == "s2_l2" or
stac_source == "aws" and collection == "sentinel-2-l2a"):
pass
else:
print(f"{collection} collection of {stac_source} not yet supported.")
# grid
if (grid != None and tile_id!= None):
if (grid == "br_states"):
br_states = load_jsons("states")
state_code = tile_id.upper()
for feature in br_states['features']:
if feature['id'] == state_code:
geom = feature['geometry']
bbox = shapely.geometry.shape(geom).bounds
geom = shapely.geometry.shape(geom["features"][0]["geometry"]) if geom["type"] == "FeatureCollection" else shapely.geometry.shape(geom)
else:
bdc_grids_data = load_jsons("grids")
selected_tile = ''
for g in bdc_grids_data['grids']:
if (g['name'] == grid):
for tile in g['features']:
if tile['properties']['tile'] == tile_id:
selected_tile = tile
geom = selected_tile['geometry']
bbox = shapely.geometry.shape(geom).bounds
geom = shapely.geometry.shape(geom["features"][0]["geometry"]) if geom["type"] == "FeatureCollection" else shapely.geometry.shape(geom)
#geometry
elif (geom):
bbox = geom.bounds
#bbox
else:
tuple_bbox = tuple(map(float, bbox.split(',')))
geom = shapely.geometry.box(*tuple_bbox)
bbox = geom.bounds
start_date = datetime.datetime.strptime(str(start_year)+'-'+str(start_month)+'-'+str(start_day), "%Y-%m-%d")
end_date = None
if all(v is not None for v in [end_year, end_month, end_day]):
end_date = datetime.datetime.strptime(str(end_year)+'-'+str(end_month)+'-'+str(end_day), "%Y-%m-%d")
elif duration_months is not None:
end_date = add_months_to_date(start_date,duration_months-1)
elif duration_days is not None and not any([end_year, end_month, end_day, duration_months]):
end_date = add_days_to_date(start_date,duration_days-1)
if end_date is None:
return print("Not provided with a valid time interval.")
periods = []
current_start_date = start_date
if duration_days:
while current_start_date <= end_date:
current_end_date = add_days_to_date(current_start_date,duration_days-1)
if current_end_date > end_date:
current_end_date = end_date
if current_start_date != current_end_date:
periods.append({
'start': current_start_date.strftime("%Y-%m-%d"),
'end': current_end_date.strftime("%Y-%m-%d")
})
current_start_date += datetime.timedelta(days=duration_days)
elif duration_months:
while current_start_date <= end_date:
current_end_date = current_start_date + dateutil.relativedelta.relativedelta(months=duration_months) - dateutil.relativedelta.relativedelta(days=1)
if current_end_date > end_date:
current_end_date = end_date
if current_start_date != current_end_date:
periods.append({
'start': current_start_date.strftime("%Y-%m-%d"),
'end': current_end_date.strftime("%Y-%m-%d")
})
current_start_date = current_start_date + dateutil.relativedelta.relativedelta(months=duration_months)
else:
periods.append({
'start': start_date.strftime("%Y-%m-%d"),
'end': end_date.strftime("%Y-%m-%d")
})
if profile=="crop_condition":
bands = ["B02","B04","B08"]
spectral_indices = ["NDVI","EVI","EVI2", "SAVI"]
elif profile=="urban_analysis":
bands = ["B02","B03","B04","B08","B11"]
spectral_indices = ["NDBI","NDVI"]
else:
spectral_indices = []
dict_collection=collection_query(
collection=collection,
start_date=start_date.strftime("%Y-%m-%d"),
end_date=end_date.strftime("%Y-%m-%d"),
bbox=bbox,
bands=bands
)
collection_name = dict_collection['collection']
if not os.path.exists(data_dir+"/"+collection):
collection_get_data(stac, dict_collection, data_dir=data_dir, stac_source=stac_source)
num_processes = multiprocessing.cpu_count()
args_for_processes = [
(period, mosaic_method, data_dir, collection_name, bands, bbox, output_dir,
duration_days, duration_months, name, geom, reference_date, projection_output, grid, tile_id, stac_source)
for period in periods
]
print(f"--- Starting parallel processing with {num_processes} processes. ---\n")
with multiprocessing.Pool(processes=num_processes) as pool:
results = pool.starmap(process_period, args_for_processes)
if(len(spectral_indices)):
calculate_spectral_indices(input_folder=output_dir,spectral_indices=spectral_indices)
if (grid_crop and not tile_id):
clip_from_grid(input_folder=output_dir, grid=grid, tile_id=tile_id)
#create_composition_json(output_dir=output_dir, collection=collection, input_scenes=scenes, ignored_scenes=[], used_scenes=[])
clean_dir(data_dir)
def process_period(period, mosaic_method, data_dir, collection_name, bands, bbox, output_dir, duration_days, duration_months, name, geom, reference_date, projection_output, grid, tile_id, stac_source):
start_date = period['start']
end_date = period['end']
process_id = os.getpid()
print(f"[Process {process_id}] Starting to process period: {start_date} to {end_date}\n")
coll_data_dir = os.path.join(data_dir+'/'+collection_name)
for i in range(0, len(bands)):
cloud_dict = get_all_cloud_configs()
cloud_list = []
band_list = []
sorted_data = []
code = stac_source.upper()+":"+collection_name
bands_cloud = [bands[i]] + [cloud_dict[code]['cloud_band']]
scenes = filter_scenes(collection_name, data_dir, geom, stac_source)
cloud = cloud_dict[code]['cloud_band']
for path in scenes:
start_dt = datetime.datetime.strptime(start_date, "%Y-%m-%d")
end_dt = datetime.datetime.strptime(end_date, "%Y-%m-%d")
filtered_files = [
f for f in os.listdir(os.path.join(coll_data_dir, path, cloud))
if (date_match := re.search(r'\d{8}', f)) and
(date_str := date_match.group()) and
len(date_str) == 8 and
date_str.isdigit() and
(file_date := datetime.datetime.strptime(date_str, "%Y%m%d")) and
start_dt <= file_date <= end_dt
]
for file in filtered_files:
date_match = re.search(r'\d{8}', file)
date_str = date_match.group()
date = datetime.datetime.strptime(date_str, "%Y%m%d")
if (reference_date):
distance_days = days_between_dates(reference_date, file.split("_")[2].split('T')[0])
pixel_count = count_pixels(os.path.join(coll_data_dir, path, cloud_dict[code]['cloud_band'], file), cloud_dict[code]['non_cloud_values'], geom)
cloud_list.append(dict(band=cloud, date=date.strftime("%Y%m%d"), distance_days=distance_days, clean_percentage=float(pixel_count['count']/pixel_count['total']), scene=path, file=''))
band_list.append(dict(band=bands[i], date=date.strftime("%Y%m%d"), distance_days=distance_days, clean_percentage=float(pixel_count['count']/pixel_count['total']), scene=path, file=''))
else:
pixel_count = count_pixels(os.path.join(coll_data_dir, path, cloud_dict[code]['cloud_band'], file), cloud_dict[code]['non_cloud_values'], geom)
if (pixel_count['total']):
cloud_list.append(dict(band=cloud, date=date.strftime("%Y%m%d"), clean_percentage=float(pixel_count['count']/pixel_count['total']), scene=path, file=''))
band_list.append(dict(band=bands[i], date=date.strftime("%Y%m%d"), clean_percentage=float(pixel_count['count']/pixel_count['total']), scene=path, file=''))
files_list = []
for path in scenes:
for band in bands_cloud:
start_dt = datetime.datetime.strptime(start_date, "%Y-%m-%d")
end_dt = datetime.datetime.strptime(end_date, "%Y-%m-%d")
filtered_files = [
f for f in os.listdir(os.path.join(coll_data_dir, path, band))
if (date_match := re.search(r'\d{8}', f)) and
(date_str := date_match.group()) and
len(date_str) == 8 and
date_str.isdigit() and
(file_date := datetime.datetime.strptime(date_str, "%Y%m%d")) and
start_dt <= file_date <= end_dt
]
for file in filtered_files:
files_list.append(dict(file=os.path.join(coll_data_dir, path, band, file)))
band_lookup, cloud_lookup = {}, {}
for f in files_list:
path = f['file']
parts = os.path.basename(path).split('_')
if (stac_source == "digitalearth-africa" and collection_name == "s2_l2a"):
date, scene, band = parts[1].split('T')[0], parts[0].lstrip('T'), parts[2].split(".")[0]
elif (stac_source == "swissdatacube" and collection_name == "s2_l2"):
date, scene, band = parts[2], parts[1], parts[3].split(".")[0]
elif (stac_source == "aws" and collection_name == "sentinel-2-l2a"):
date, scene, band = parts[2], parts[1], map_band_name(parts[3].split(".")[0])
elif (stac_source == "bdc" and collection_name == "S2_L1C_BUNDLE-1" or
stac_source == "bdc" and collection_name == "S2_L2A-1"):
date, scene, band = parts[2].split('T')[0], parts[5].lstrip('T'), parts[1]
#else:
# date, scene, band = parts[1].split('T')[0], parts[0].lstrip('T'), parts[2].split(".")[0]
cloud_lookup[(date, scene)] = path
band_lookup[(band, date, scene)] = path
for item in band_list:
item['file'] = band_lookup.get((item['band'], item['date'], item['scene']), '')
for item in cloud_list:
item['file'] = cloud_lookup.get((item['date'], item['scene']), '')
if (mosaic_method=='lcf'):
sorted_data = sorted(band_list, key=lambda x: x['clean_percentage'], reverse=True)
cloud_sorted_data = sorted(cloud_list, key=lambda x: x['clean_percentage'], reverse=True)
if (mosaic_method=='chrono'):
sorted_data = sorted(band_list, key=lambda x: x['date'])
cloud_sorted_data = sorted(cloud_list, key=lambda x: x['date'])
if (mosaic_method=='ctd'):
sorted_data = sorted(band_list, key=lambda x: x['distance_days'])
cloud_sorted_data = sorted(cloud_list, key=lambda x: x['distance_days'])
reproject_data = reproject_tifs(sorted_data=sorted_data, cloud_sorted_data=cloud_sorted_data, data_dir=data_dir, projection_output=projection_output)
sorted_data = reproject_data['reprojected_images']
cloud_sorted_data = reproject_data['reprojected_cloud_images']
if (i==0):
ordered_lists = merge_scene_provenance_cloud(sorted_data, cloud_sorted_data, scenes, collection_name, bands[i], data_dir, stac_source, start_date, end_date)
else:
ordered_lists = merge_scene(sorted_data, cloud_sorted_data, scenes, collection_name, bands[i], data_dir, stac_source, start_date, end_date)
filename = sorted_data[0]['file'].split('/')[-1]
if (collection_name =='S2_L2A-1'):
baseline_number = filename.split("_N")[1][0:4]
else:
baseline_number = 0
band = bands[i]
if not os.path.exists(output_dir):
os.makedirs(output_dir)
collection_prefix = collection_name.split("-")[0].upper()
name_upper = name.upper()
date_range = f"{str(start_date).replace('-', '')}_{str(end_date).replace('-', '')}"
current_band = bands[i]
if duration_months:
duration_str = f"-{duration_months}M"
elif duration_days:
duration_str = f"-{duration_days}D"
else:
duration_str = ""
file_name = f"{collection_prefix}-{name_upper}{duration_str}-{current_band}_{date_range}"
cloud_file_name = f"{collection_prefix}-{name_upper}{duration_str}_{cloud}_{date_range}"
provenance_file_name = f"{collection_prefix}-{name_upper}{duration_str}-PROVENANCE_{date_range}"
output_file = os.path.join(output_dir, f"raw-{file_name}.tif")
if i == 0:
cloud_data_output_file = os.path.join(output_dir, f"cloud_data_raw-{file_name}.tif")
provenance_output_file = os.path.join(output_dir, f"provenance_raw-{file_name}.tif")
datasets = [rasterio.open(file) for file in ordered_lists['merge_files']]
extents = get_dataset_extents(datasets)
merge_tifs(tif_files=ordered_lists['merge_files'], output_path=output_file, band=band, path_row=name, extent=extents)
if (i==0):
merge_tifs(tif_files=ordered_lists['provenance_merge_files'], output_path=provenance_output_file, band=band, path_row=name, extent=extents)
merge_tifs(tif_files=ordered_lists['cloud_merge_files'], output_path=cloud_data_output_file, band=cloud_dict[code]["cloud_band"], path_row=name, extent=extents)
clip_raster(input_raster_path=output_file, output_folder=output_dir, clip_geometry=geom, projection_output=projection_output, output_filename=file_name+".tif", grid=grid, tile_id=tile_id)
if (i==0):
clip_raster(input_raster_path=cloud_data_output_file, output_folder=output_dir, clip_geometry=geom,projection_output=projection_output, output_filename=cloud_file_name+".tif", grid=grid, tile_id=tile_id)
clip_raster(input_raster_path=provenance_output_file, output_folder=output_dir, clip_geometry=geom, projection_output=projection_output, output_filename=provenance_file_name+".tif", grid=grid, tile_id=tile_id)
fix_baseline_number(input_folder=output_dir, input_filename=file_name, baseline_number=baseline_number)
generate_cog(input_folder=output_dir, input_filename=file_name, compress='DEFLATE')
if (i==0):
generate_cog(input_folder=output_dir, input_filename=cloud_file_name, compress='DEFLATE')
generate_cog(input_folder=output_dir, input_filename=provenance_file_name, compress='DEFLATE')
clean_dir(data_dir=data_dir,date_interval=str("-"+str(start_date).replace("-", "")+'_'+str(end_date).replace("-", "")))