[GIS] Creating composite band mosaic of sentinel 2 data with arcpy

arcpycompositeiteratorremote sensingsentinel-2

I previously useda simple script that went through many Landsat satellite images, and for each folder it found it made a composite band raster. The script was as follows:

import arcpy
import os

arcpy.env.workspace = (folder containing lots of Landsat downloads, all of which in a separate folder)
outws = (output folder)
folders = arcpy.ListWorkspaces()

for folder in folders:
    arcpy.env.workspace = folder
rasters = arcpy.ListRasters("*.tif")
name = os.path.join(outws, rasters[1].split("_")[0] + ".tif")
arcpy.CompositeBands_management(rasters, join)

print "Processing complete"

This was very useful as it saves lots of time. I now want to adapt it to work with the new Sentinel 2 satellite images. There are some major differences though. For Landsat, there was a folder which contained 8 bands, all of which covered the same area. For Sentinel 2 the bands are split up into different subsets (the image I have here is 14 subsets) so somehow I would like to make the script navigate to where the subsets are stored, this is a folder called GRANULE. Within this folder there is a folder for each subset, the bands that I would like to combine are stored in a subfolder called IMG_DATA.

So, I would like the script to find the 14 folders called IMG_DATA, make a composite band for each folder, and finally mosaic them into one raster. Does someone know how I could adapt my script to use this?

Best Answer

you can replace the list raster with glob.glob in order to get the list of your raster : 

import glob, arcpy

list_composites = []
list_images = glob.glob("path_to_your_image.SAFE\GRANULE\*\IMG_DATA\*B01.jp2") #first band in each folder
for image in list_images:
    rasters = glob.glob(image[:-6]+"*") #all bands in one folder
    arcpy.CompositeBands_management(rasters[1:], image[:-6]+"join.tif")#I removed the first band (aerosol, 60m). Note that composite bands takes the pixel size of the first band in the list.
    list_composites.append(image[:-6]+"join.tif")

then you can use your list_composites to create a mosaic or a raster catalog.

Note that i) with SNAP, you can create the composite and extract a subset without unzipping, ii) with gdalbuildvrt you can create a virtual file to avoid duplicating your data (Sentinel-2 is larger than Landsat-8, this makes big files).

Related Question