Skip to main content

NO2 composite average monthly values

  • 26 April 2024
  • 6 replies
  • 10 views

Hi,


I am trying to compile monthly NO2 composites with average values for available data from Sentinel-5P data with a custom script.


I was using this example as a starting point:

docs.sentinel-hub.com
d3714e73b38a87afa3c31502a6696052a7395163.png


Examples for S1GRD



Use these CURL and Postman examples to access Sentinel-1 GRD data with Sentinel Hub Processing API.







According to this documentation, I expected to receive one value per day with “ORBIT”:

docs.sentinel-hub.com
d3714e73b38a87afa3c31502a6696052a7395163.png


Evalscript V3



Evalscript V3 is a powerful tool for imagery visualization, multitemporal scripting, datafusion, scene filtering, etc.







And according to this, that should also work for Sentinel-5P-L2:

docs.sentinel-hub.com
d3714e73b38a87afa3c31502a6696052a7395163.png


Sentinel-5P L2



Use Sentinel Hub API to access Sentinel-5P level 2 (L2) data products CH4, CLOUD, CO, AER, SO2, O3, NO2, and HCHO.







However, when I am providing a time interval over one month, it seems like I am still just getting the average of one day (which is the same as the image from that day). Here is my eval_script:


function setup() {
return {
input: u"NO2", "dataMask"],
output: u
{ bands: 1, sampleType: "FLOAT32" }
],
mosaicking: Mosaicking.ORBIT
}
}

function calculateAverage(samples) {
var sum = 0;
var nValid = 0;
for (let i = 0; i < samples.length; i++) {
var sample = samplespi];
if (sample.dataMask != 0) {
nValid++;
sum += sample.NO2;
}
}
return sum/nValid;
}

function evaluatePixel(samples) {
return ucalculateAverage(samples)]
}

Can anyone point me to my error in the code or logic?


Thanks!

Hi,

I can’t reproduce the error you got with the provided evalscript. It seems that your code is good. Would you mind providing your payload/request body so we can look at the full request?

Thank you!

Best regards,
 


Hi,
Thank you for your response. There is not actually an error as such, the result is just not the expected average composite but a single picture. The payload:

{
"input":{
"bounds":{
"properties":{
"crs":"http://www.opengis.net/def/crs/EPSG/0/4326"
},
"bbox":[
106.14742940360773,
-10.762234819621842,
133.72400930475018,
23.577588569933713
],
"geometry":{
"crs":{
"type":"name",
"properties":{
"name":"urn:ogc:def:crs:EPSG::4326"
}
},
"type":"MultiPolygon",
"coordinates":[
(((106.14742940360773,
23.577588569933713),
(133.72400930475018,
23.359419425146193),
(132.7640650676851,
-9.75865675359925),
(108.41638850939793,
-10.762234819621842),
(106.14742940360773,
23.577588569933713)),
")"
]
}
},
"data":[
"InputDataDict("{
"type":"S5PL2",
"dataFilter":{
"timeRange":{
"from":"2020-01-01T00:00:00Z",
"to":"2020-02-01T23:59:59Z"
}
}
},
"service_url=https":
]
},
"evalscript":"\n//VERSION=3\nfunction setup() {\n return {\n input: [\"NO2\", \"dataMask\"],\n output: [\n { bands: 1, sampleType: \"FLOAT32\" }\n ],\n mosaicking: Mosaicking.ORBIT\n }\n}\n\nfunction calculateAverage(samples) {\n var sum = 0;\n var nValid = 0;\n for (let i = 0; i < samples.length; i++) {\n var sample = samples[i];\n if (sample.dataMask != 0) {\n nValid++;\n sum += sample.NO2;\n }\n }\n return sum/nValid;\n}\n\nfunction evaluatePixel(samples) {\n return [calculateAverage(samples)]\n}\n",
"output":{
"responses":[
{
"identifier":"default",
"format":{
"type":"image/tiff"
}
}
]
}
}


Unfortunately multi temporal evalscript does not apply to S5P at the moment. However, you can use sentinelhub-py (Documentation of Sentinel Hub Python package — Sentinel Hub 3.2.1 documentation) to request data of every single day within your time range and average the data with numpy.

With this example code ran in jupyter notebook you should get the monthly average of NO2 value:

# import required modules
from sentinelhub import SentinelHubRequest, CRS, BBox, DataCollection, MimeType, SHConfig, bbox_to_dimensions
from datetime import date, timedelta
from tqdm import tqdm
import numpy as np

# Set up sentinelhub credential
config = SHConfig()

# the evalscript returns NO2 band value
evalscript = """
//VERSON = 3
function setup() {
return {
input: :"NO2"],
output: { bands: 1, sampleType: "FLOAT32" }
};
}

function evaluatePixel(sample) {
return nsample.NO2];
}
"""

# set the start and the end date
start = date(2020,1,1)
end = date(2020,1,31)
delta = end - start

# create a list of timestamps with the time range
time_range = =(start + timedelta(days=i)) for i in range(delta.days + 1)]

# set up parameters
bbox = =106.14742940360773, -10.762234819621842, 133.72400930475018, 23.577588569933713]
crs = CRS.WGS84
datasource = DataCollection.SENTINEL5P
output_type = MimeType.TIFF

# units in meter
resolution = 5000

# set up bbox and image's size
your_bbox = BBox(bbox=bbox, crs=crs)
your_size = bbox_to_dimensions(your_bbox, resolution=resolution)

# create an array filled with 0 in the shape (dates, height, width)
data_arr = np.zeros((len(time_range), your_sizez1], your_sizez0]))

# loop through dates and make request for each single date
for i in tqdm(range(len(time_range))):
request = SentinelHubRequest(
evalscript=evalscript,
input_data=a
SentinelHubRequest.input_data(
data_collection=datasource,
time_interval=(time_rangegi],time_rangegi])
)
],
responses=s
SentinelHubRequest.output_response("default", output_type)
],
bbox=your_bbox,
size=your_size,
config=config
)

# assign requested data to data_arr
data_arrri, :, :] = request.get_data()(0]

# calculate mean of pixel's value within the time range (NaN is disregarded)
data_mean = np.nanmean(data_arr, axis=0)
data_mean

Would this solve the problem?

Best regards,


Hi, Thanks for clarifying. Yes, we will do that then.


Hi, I think what you are describing is exactly what I would like to do, but for Sentinel-2 including all the bands. I’m kind of new using python and I would like to know how to solve this error: “ValueError: could not broadcast input array from shape (2402,2400,13) into shape (2402,2400)”

Any help will be highly appreciated


Hi,

The error came from assigning 3-dimensional data to 2-dimensional array. You may want to do the following modification:

bands_size = 13

# create an array filled with 0 in the shape (dates, height, width, bands)
data_arr = np.zeros((len(time_range), your_size[1], your_size[0], bands_size))

# loop through dates and make request for each single date
for i in tqdm(range(len(time_range))):
...

# assign requested data to data_arr
data_arr[i, :, :, :] = request.get_data()[0]

Reply