Zonal Statistics - obtain min and max value for the whole year for each land type

Microwave Toolbox
#1

I would like to obtain minimum and maximum values of the Sentinel-1 backscatter coefficient for each type of the land usage and export them to CSV/text file.

Sample output from my analysis will looks as follows:

| Type | Min* | Max* |
| Forests | Val1 | Val2 |
| Meadows | Val3 | Val3 |

where * - min and max values for year 2015.

I thought about preparing shapefile which will define extent of each type (forests, meadows) and then calculate minimum and maximum value for each area. However, I’m not sure if there’s any counterpart of QGIS Zonal Statistics plugin in ESA SNAP Desktop?
(I’ve checked topic Compute zonal statistics from Sentinel-2 VI’s , but I can’t find answer to my question)

#2

I have an idea, I didn’t implement, but might be it is helpful,
Convert the result (S1) to geotiff, and then open the result in ArcMap, create shfile polygons and extract the values of any land type,

#3

Thank you, but I would like to avoid this step. Since Sentinel-1 products are quite large and GeoTIFF is not optimal format to store them. As I mentioned in my question, there’s appropriate tool in QGIS, I don’t need to use ArcGIS (which is Windows only tool). And ideally, I would like to stay in ESA SNAP Desktop/Python scripts.

#4

The max/min/etc. values can be achieved quite easily with the BandMaths operator, which you can iterate over the timeseries.

GPT iteration comparison between a new band(band_00) and old max/min value(band_min/band_max):

I am using BandMerge as a source, where I initially load two separate BEAM-DIMAP files, one with the stats, the other with input.
BandMerge the BandMaths operators at the end of GPT to get single file or write multiple files for each.

<node id="BandMaths_max">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="BandMerge"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
        <targetBand>
          <name>sigma0_max</name>
          <type>float32</type>
          <expression>${band_00} &gt;= ${band_max} ? ${band_00} : ${band_max}</expression>
          <description/>
          <noDataValue>0.0</noDataValue>
        </targetBand>
      </targetBands>
      <variables/>
    </parameters>
  </node>

  <node id="BandMaths_min">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="BandMerge"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
        <targetBand>
          <name>sigma0_min</name>
          <type>float32</type>
          <expression> ${band_00} &lt;= 0.0005 ? ${band_min} : ${band_min} &lt;= 0.0005 ? ${band_00} : (${band_min} &lt;= ${band_00} ? ${band_min} : ${band_00}) </expression>
          <description/>         
          <noDataValue>0.0</noDataValue>
        </targetBand>
      </targetBands>
      <variables/>
    </parameters>
  </node>

For the classification i would recommend similar approach, just define an interval for corresponding values, sadly I can’t be helpful with that one.

2 Likes
#5

Thank you.
I understand (or hope so) your general idea, but I have problem with using provided by you .xml.

I thought that if I go to Tools --> Graph Builder —> Load xml and select your xml, I will see new blocks, but I wasn’t.
Could you be so kind and suggest any tutorials or provide me some hints of how can I use your .xml or which tutorial should I study?

Moreover, I don’t need to do classification, sorry if it wasn’t clear.
I will have a shapefile

#6

The Graph is just 2 nodes (boxes) extracted, with all the variables used its probably better to use the gpt from the console, you can call it from a script java, python as well. There you can define the parameters as you are only using one band from the source file at a time so loop is a nice thing to do there.

This Graph should be complete now to give you an initial overview:
graph_stats_simple.xml (3.3 KB)

For the classification you mentioned the forest and meadows so I was aiming at that, if that is not the aim, the simpler.

With the shapefile it makes sense, I believe raster conversion via QGIS or GDAL will do nicely.

Good luck, hope this helps a bit.

Petr

#7

I just wonder if it is possible to do that without scripting / using GUI tools only. I have some experience with Python, I will try to it through Python.

Could you tell me if you suggest to convert vector file (.shp) to raster format and use it in Raster Calculator / Band Math in SNAP?
Thank you in advance for response.

#8

Vector file to raster can be done easily with QGIS, rasterize is the function I believe.

Otherwise you can load shapefiles directly into the SNAP Toolbox, then you can convert/handle them as masks. Even better as the shapefiles can be integrated into the BEAM-DIMAP then you can use it in BandMaths as you wanted.

These are just a few possible options.

#9

Hi Suprd,
In the xml file, I didn’t understood the how do you write the bandmath expression
${band_00} >= ${band_max} ? ${band_00} : ${band_max} how does this generate the band maximum, where is the looping here?

${band_00} <= 0.0005 ? ${band_min} : ${band_min} <= 0.0005 ? ${band_00} : (${band_min} <= ${band_00} ? ${band_min} : ${band_00}) how does this generate the band minimum

For suppose I want to generate the 70th percentile value how to do it?

#10

I tried to use your xml for my raster file. It is giving the following error.

image.png1087×30 2.15 KB