NO2 composite average monthly values

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 [sample.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_size[1], your_size[0]))

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

    # assign requested data to data_arr
    data_arr[i, :, :] = 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]