Thursday 11 January 2018

R script for generating multi-line graphs from .csv files

Recently I got reviewer's comments on a paper and they asked me to remake my graphs using R. At the beginning I got a bit annoyed that I had to redo them, but the results looks great!! :)  Therefore, I thought about documenting the code within my blog for future reference.

Let assume that you have the following .csv file and you would like to create a multi-line graph (please note that Excel also allows to export spreadsheet into .csv file format).

Content of a file table4.csv
The following R script creates a multi-line graphs. The input variables at the beginning can be adjusted accordingly.

### User defined Variables ###
# The title that appears on the top of the graph
mainTitle  = "Recall"
# The tile of the x-axis
xAxisTitle = "Distance (m)"
# The title of the y-axis
yAxisTitle = "Recall"
# The name of the imported .csv file
fileName = "table4.csv"
# The name of the image where the graph will be stored
exportFileName = "table4.png"

# Read the .csv file
data <- read.csv(file=fileName,row.names=1, head=TRUE,sep=",",check.names=FALSE)

# extracts two lists with the labels of the axes
rowNames = row.names(data)
colNames = colnames(data)

png(filename=exportFileName)

# plot graph with the given titles
matplot(t(data), type="l", lty=1, lwd=2, main= mainTitle, xlab=xAxisTitle, ylab=yAxisTitle)

# modify numbering\labels of the axes to agree with the inserted table of the .csv file
axis(1, at=1:length(colNames), lab=colNames)

# Add legend to the table
legend("topleft", inset=0.01, legend=rowNames, col=c(1:6), bg= ("white"), horiz=F, lty=1)
dev.off()

The result is the following graph:
The graph created using the R script

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.