I want to create composite band image with sentinel 2 data through arcpy.

396
5
01-19-2021 09:26 AM
SiddharthaDhara
New Contributor II

As we know sentinel satellite data have many sub folder. I want to create composite band of each datasets with the band 2, band 3, band 4, band 8 which are located in one folder(10m) and band 11, band 12 which are located in other folder(20m).I am making a arcpy code for that but code was not run correctly. Here is the code

 

 

import arcpy, os

#arcpy.env.workspace = "L:\Arcpy\Data\Satellite_image_raw\New folder"

basepath = "L:\Arcpy\Data\Satellite_image_raw\New folder"

endswith = ("_B02_10m.jp2", "_B03_10m.jp2", "_B04_10m.jp2", "_B08_10m.jp2", "_B11_20m.jp2", "_B12_20m.jp2")

output = "L:\Arcpy\Data\Layer_stack"

rasterlist = []

#folders = arcpy.ListWorkspaces()

#for folder in folders:
#arcpy.env.workspace = folder
for root, dirs, files in os.walk(basepath):
     for rasterlayer in files:
           if rasterlayer.endswith(endswith):
              #print(os.path.join(root, rasterlayer))
              rasterlist.append(rasterlayer)
              print "image:" + str(rasterlist)
name = os.path.join(output, rasterlist[1].split("_")[0] + ".img")
arcpy.CompositeBands_management(rasterlist, name)
print(rasterlist)

 

0 Kudos
5 Replies
DanPatterson
MVP Esteemed Contributor

This will get you started... you need to raw encode your paths. (the little "r" )

p = "L:\Arcpy\Data\Satellite_image_raw\New folder"  # --- bad
  File "<ipython-input-3-74ab7250855b>", line 1
    p = "L:\Arcpy\Data\Satellite_image_raw\New folder"
       ^
SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 33-34: malformed \N character escape


p = r"L:\Arcpy\Data\Satellite_image_raw\New folder"  # --- good

... sort of retired...
0 Kudos
SiddharthaDhara
New Contributor II

import glob, arcpy, os

arcpy.env.workspace = r"D:\Composite_band\New folder"

folders = arcpy.ListWorkspaces()

list_images = glob.glob (r"D:\Composite_band\New folder\S2A_MSIL2A_20210104T044201_N0214_R033_T45QVD_20210104T063322.SAFE\GRANULE\L2A_T45QVD_A028916_20210104T044202\IMG_DATA\R10m\T45QVD_20210104T044201_AOT_10m.jp2") #first 10m folder in each folder
list_images1 = glob.glob (r"D:\Composite_band\New folder\S2A_MSIL2A_20210104T044201_N0214_R033_T45QVD_20210104T063322.SAFE\GRANULE\L2A_T45QVD_A028916_20210104T044202\IMG_DATA\R20m\T45QVD_20210104T044201_AOT_20m.jp2") #Second 20m folder in each folder

output = r"D:\Composite_band\New folder\Layer_stack"

for folder in folders:
      arcpy.env.workspace = folder
      for image in list_images:
            rasters = glob.glob(image[1:-5]) #Selected bands in 10m folder
            for image1 in list_images1:
                  rasters1 = glob.glob(image1[8:10]) #Selected bands in 20m folder
                  rasters.append(rasters1)
                  name = os.path.join(output, rasters[0].split("_")[0] + ".img")
                  arcpy.CompositeBands_management(rasters,name)
                  print(rasters)

I was trying as you said but output is not created. It is gave this error Traceback (most recent call last):
File "D:\Composite_band\composite_band_script_1.py", line 32, in <module>
name = os.path.join(output, rasters[0].split("_")[0] + ".img")
AttributeError: 'list' object has no attribute 'split'

0 Kudos
DanPatterson
MVP Esteemed Contributor

you had better throw some print statements in there because you are assuming rasters[0] is indeed a raster and not something else and you won't be able to get "name" until you are sure it is a single file name and not a raster, which I suspect it is


... sort of retired...
0 Kudos
SiddharthaDhara
New Contributor II

This time I am trying this way 

import arcpy, os

arcpy.env.workspace = ("D:\Composite_band\New folder")

basepath = ("D:\Composite_band\New folder")

endswith = ("_B02_10m.jp2", "_B03_10m.jp2", "_B04_10m.jp2", "_B08_10m.jp2", "_B11_20m.jp2", "_B12_20m.jp2")

output = ("D:\Composite_band\Layer_stack")

folders = arcpy.ListWorkspaces()

for folder in folders:
      arcpy.env.workspace = folder
      for root, dirs, files in os.walk(basepath):
            for file in files:
                  if file.endswith(endswith):
                  print(os.path.join(root, file))
                  image = []
                  image = file
                  print(image)
                  name = os.path.join(output, image[1].split("_")[0] + ".img")
                  arcpy.CompositeBands_management(image, name)
                  print(image)

But this time code is able to print first selected band of the first folder after that it gives error like this:  D:\Composite_band\New folder\S2A_MSIL2A_20210104T044201_N0214_R033_T45QVD_20210104T063322.SAFE\GRANULE\L2A_T45QVD_A028916_20210104T044202\IMG_DATA\R10m\T45QVD_20210104T044201_B02_10m.jp2
T45QVD_20210104T044201_B02_10m.jp2
Traceback (most recent call last):
File "D:\Composite_band\composite_band_script_1.py", line 34, in <module>
arcpy.CompositeBands_management(image, name)
File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 13595, in CompositeBands
raise e
ExecuteError: ERROR 000271: Cannot open the input datasets
Failed to execute (CompositeBands).

It means after printing first band the code is not able to open the input folder.

0 Kudos
DanPatterson
MVP Esteemed Contributor

you changed the code

     for root, dirs, files in os.walk(basepath):

I suspect you are trying to piece something together without fully understanding the implications.

I would suggest using the builtin tools to see if you can get a workflow established before continuing down the path of trying to script it.

Good luck


... sort of retired...
0 Kudos