Thank you @chung.horng for the prompt reply. Here is my code converted from a Jupyter Notebook. If there is a more appropriate way for me to upload the code, please let me know for future issues.
get_ipython().run_line_magic(‘matplotlib’, ‘inline’)
import datetime
import matplotlib.pyplot as plt
from aenum import MultiValueEnum
from matplotlib.colors import BoundaryNorm, ListedColormap
import geopandas as gpd
import math
from sentinelhub import CRS, BBox, DataCollection, SHConfig
from eolearn.core import EOWorkflow, FeatureType, LoadTask, OutputTask, SaveTask, linearly_connect_tasks, MergeFeatureTask
from eolearn.io import SentinelHubDemTask, SentinelHubEvalscriptTask, SentinelHubInputTask, ExportToTiffTask
In[2]:
In case you put the credentials into the configuration file you can leave this unchanged
CLIENT_ID = “”
CLIENT_SECRET = “”
In[3]:
config = SHConfig()
if CLIENT_ID and CLIENT_SECRET:
config.sh_client_id = CLIENT_ID
config.sh_client_secret = CLIENT_SECRET
if config.sh_client_id == “” or config.sh_client_secret == “” or config.instance_id == “”:
print(“Warning! To use Sentinel Hub services, please provide the credentials (client ID and client secret).”)
In[4]:
shapefile = “/home/nedhorning/Abe/LCL_Farm_Parcels/testParcelsUTM18N.shp”
eopatch_path = “/home/nedhorning/Abe/sentinel_hub_data/”
#eopatch_dir=“eopatch”
output_allbands_filename = “testUTM21S.tif”
output_ndvi_filename = “NDVI_UTM21S.tif”
epsg = “EPSG:32618”
time interval of downloaded data
time_interval = (“2022-07-01”, “2022-08-30”)
resolution of the request (in metres)
resolution = 30
expansion_dist = 210 # Distance in meters to expance each side of the bounding rectangle
maximal cloud coverage (based on Sentinel-2 provided tile metadata)
#maxcc = 0.1
time difference parameter (minimum allowed time difference; if two observations are closer than this,
they will be mosaicked into one observation)
time_difference = datetime.timedelta(hours=2)
In[5]:
gdf = gpd.read_file(shapefile)
In[6]:
bands_evalscript = “”"
//VERSION=3
function setup() {
return {
input: [{
bands: ["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B10", "BQA", "QA_RADSAT", "SR_QA_AEROSOL", "ST_QA", "ST_TRAD", "ST_URAD", "ST_DRAD", "ST_ATRAN", "ST_EMIS", "ST_EMSD", "ST_CDIST"]
}],
output: {
id: "ms_data",
bands: 19,
sampleType: SampleType.UINT16
}
}
}
function evaluatePixel(sample) {
return [10000 * sample.B01, 10000 * sample.B02, 10000 * sample.B03, 10000 * sample.B04, 10000 * sample.B05, 10000 * sample.B06, 10000 * sample.B07, sample.B10, sample.BQA, sample.QA_RADSAT, sample.SR_QA_AEROSOL, sample.ST_QA, sample.ST_TRAD, sample.ST_URAD, sample.ST_DRAD, sample.ST_ATRAN, sample.ST_EMIS, sample.ST_EMSD, sample.ST_CDIST]
}
“”"
In[7]:
indices_evalscript = “”"
//VERSION=3
function setup() {
return {
input: ["B04","B05"],
output:[{
id: "indices",
bands: 1,
sampleType: SampleType.INT16
}]
}
}
function evaluatePixel(sample) {
let ndvi = Math.round(index(sample.B05, sample.B04) * 10000);
return {
indices: [ndvi]
};
}
“”"
In[8]:
this will add two indices: ndvi, gci
add_indices = SentinelHubEvalscriptTask(
features=[(FeatureType.DATA, “indices”)],
evalscript=indices_evalscript,
data_collection=DataCollection.LANDSAT_OT_L2,
resolution=resolution,
#maxcc=maxcc,
time_difference=time_difference,
config=config,
max_threads=3,
)
In[9]:
this will add all Sentinel bands as DN (int) not reflectance (float)
add_ms_bands = SentinelHubEvalscriptTask(
features=[(FeatureType.DATA, “ms_data”)],
evalscript=bands_evalscript,
data_collection=DataCollection.LANDSAT_OT_L2,
resolution=resolution,
#maxcc=maxcc,
time_difference=time_difference,
config=config,
max_threads=3,
)
In[10]:
ecopatch_dir_list = []
for i in range(len(gdf)):
print(f"\rCalculating parcel {i+1} of {len(gdf)}", end=‘’, flush=True)
# Get the current polygon
polygon = gdf.loc[i, “geometry”]
# Extract the polygon's coordinates
x_coords, y_coords = polygon.exterior.coords.xy
# Get the filename by concatenating the "farm" and "name" attributes
eopatch_dir = gdf.loc[i, 'FARM'] + "_" + gdf.loc[i, 'Name']
ecopatch_dir_list.append(eopatch_dir)
xmin = int(min(x_coords) - expansion_dist)
xmax = int(max(x_coords) + expansion_dist)
ymin = int(min(y_coords) - expansion_dist)
ymax = int(max(y_coords) + expansion_dist)
# Adjust BBox so it is evenly divisible by resolution - if not then final resolution is constrained by BBox
if (((ymax - ymin) % resolution) != 0):
ymax = ymax - ((ymax - ymin) % resolution) + resolution
if (((xmax - xmin) % resolution) != 0):
xmax = xmax - ((xmax - xmin) % resolution) + resolution
# region of interest
roi_coordinates = [xmin, ymin, xmax, ymax]
roi_bbox = BBox(roi_coordinates, crs=CRS(epsg))
save = SaveTask(eopatch_path + "numpyLandsat/", overwrite_permission=2, compress_level=0)
output_task = OutputTask(eopatch_dir)
workflow_nodes = linearly_connect_tasks(add_indices, add_ms_bands, save, output_task)
workflow = EOWorkflow(workflow_nodes)
result = workflow.execute(
{
workflow_nodes[0]: {"bbox": roi_bbox, "time_interval": time_interval},
workflow_nodes[-2]: {"eopatch_folder": eopatch_dir},
}
)
