Friday 5 January 2018

Extracting, geocorrecting and exporting a band into GeoTiff; Sentinel 3 SLSTR instrument

I have recently started working with Sentinel 3 data and even though they are really cool data, sometimes it is a hassle to find the appropriate tutorial and\or the information you need. It took me a while to automated the process of extracting a band and saving it into .tif. Therefore, I thought more that people will be interested into the solution.

For this tutorial you will need Python and SNAP installed.


You may download images from the SLSTR instrument using the following link:
https://coda.eumetsat.int/
Please note that the following approach does not work for the OLCI instrument. SNAP prints an error message. I hope it will be fixed soon.

Also do not try gdal_translate, it gives you an output and you may think it is ok at the beginning but looking deeper into it, the coordinate locations of the pixels are not correct. And this happens because the images are not georectificated. The coordinates of each pixel are stored into the lat  and lot. The best approach then is to use SNAP whose purpose is to manage Sentinel data. 


You can extract and reproject\georectificate a band using the Graph Builder (found in the Tools menu) as shown in the following image:

Figure 1: Graph Builder of SNAP

At first you need to select a band to avoid geocorreting (reprojecting) all the bands, because it is a time consuming tasks. The output of the BandSelect node is a 3 band image, because the lon and lat bands are preserved and are necessary components for the reprojection. The Reprojection command georectificate the image. But all three bands are preserved into its output. So we need to select again the band of our interest in order to drop the lat and lon bands.


If you want to automate the process, you can export the graph into an .xml file and run it into a terminal using the "gpt" command. You may modify the .xml file and add variables using the following: ${variableName}. Then you can define it from the terminal this way: -PvariableName="variableName". An example follows.

The following .xml files was exported from the graph illustrated in Figure 1 and modified to add the variables "in", "out" and "band".


<graph id="Graph">
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <file>${in}</file>
    </parameters>
  </node>
  <node id="BandSelect">
    <operator>BandSelect</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <selectedPolarisations/>
      <sourceBands/>
      <bandNamePattern>${band}</bandNamePattern>
    </parameters>
  </node>
  <node id="Reproject">
    <operator>Reproject</operator>
    <sources>
      <sourceProduct refid="BandSelect"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <wktFile/>
      <crs>GEOGCS["WGS84(DD)", 
  DATUM["WGS84", 
    SPHEROID["WGS84", 6378137.0, 298.257223563]], 
  PRIMEM["Greenwich", 0.0], 
  UNIT["degree", 0.017453292519943295], 
  AXIS["Geodetic longitude", EAST], 
  AXIS["Geodetic latitude", NORTH]]</crs>
      <resampling>Nearest</resampling>
      <referencePixelX/>
      <referencePixelY/>
      <easting/>
      <northing/>
      <orientation/>
      <pixelSizeX/>
      <pixelSizeY/>
      <width/>
      <height/>
      <tileSizeX/>
      <tileSizeY/>
      <orthorectify>false</orthorectify>
      <elevationModelName/>
      <noDataValue>NaN</noDataValue>
      <includeTiePointGrids>true</includeTiePointGrids>
      <addDeltaBands>false</addDeltaBands>
    </parameters>
  </node>
  <node id="BandSelect(2)">
    <operator>BandSelect</operator>
    <sources>
      <sourceProduct refid="Reproject"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <selectedPolarisations/>
      <sourceBands>${band}</sourceBands>
      <bandNamePattern>${band}</bandNamePattern>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="BandSelect(2)"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <file>${out}</file>
      <formatName>GeoTIFF</formatName>
    </parameters>
  </node>
  <applicationData id="Presentation">
    <Description/>
    <node id="Read">
            <displayPosition x="37.0" y="134.0"/>
    </node>
    <node id="BandSelect">
      <displayPosition x="163.0" y="141.0"/>
    </node>
    <node id="Reproject">
      <displayPosition x="309.0" y="138.0"/>
    </node>
    <node id="BandSelect(2)">
      <displayPosition x="451.0" y="154.0"/>
    </node>
    <node id="Write">
      <displayPosition x="611.0" y="180.0"/>
    </node>
  </applicationData>
</graph>

You may run the script as follow:

       
 gpt ExtractReprojectBandfromS3_SLSTR.xml -Pin=dir\Sen3SLSTRfile.nc +  -Pband="band_name" -Pout="out.tiff"

Acknowledgements
This tutorial was funded under the SEO-DWARF project of H2020 RISE
Written during my secondment as researcher from Cyprus University of Technology to Planetek Hellas.

No comments:

Post a Comment