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.

Wednesday, 9 August 2017

The structure of the LAS1.3 file format used to store full-waveform LiDAR data

There are a few LiDAR file formats but the LAS1.3 was the first format to contain FW data and it is the file format supported by DASOS (http://miltomiltiadou.blogspot.gr/2015/03/las13vis.html). According to the LAS1.3 file specifications, a .LAS file contains information about both discrete and FW LiDAR data, with the waveform packets attached to discrete returns and saved either internally at the end of the .LAS file or externally in a .WVS file.

As shown at the following Figure, the .LAS file is divided into four sections.


 A brief explanation of each section is given here:
  • The Header contains general information about the entire flightline. For example, it includes the maximum scan angle used during the flight, whether the waveform packets are recorded internally or externally  and the number of Variable Length Records (VLR).
  •  Regarding the VLR, which contain arbitrary "extension" data blocks, the most important information given is the waveform packet descriptors that contain essential information on how to read the waveform packets (i.e. an ID, the number of wave samples and the size of each intensity in bits).
  •  The Point Data Records are the discrete points and the waveforms are associated with first return discrete points. Each Point Data Record has a spatial location, an intensity and optionally a pointer to a waveform packet as well as the ID of the corresponding waveform packet descriptor.
  •  The waveform packets is a list of intensities and they are either saved internally into the Extended Variable Length Records section of the .LAS file or inside an external .WVS file. Starting from the associated first return point, the spatial locations of the waveform packet (wave sample intensity) are calculated by adding an offset defined in the associated Point Data Record.
The full LAS1.3 file specifications are available here: http://www.asprs.org/a/society/committees/standards/LAS_1_3_r11.pdf

Please note that this text was taken from the EngD thesis of Milto Miltiadou, which was submitted to University of Bath in 2017.
   

Wednesday, 4 January 2017

Statistics in C++ (Mean, Median and Standard Deviation)

Just a few simple statistics. I found this function extremely useful for my software DASOS. I used it for calculating features that may characterise trees. :)

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <cmath>

typedef struct Statistics {
   double mean;
   double median;
   double stdev;
} Statistics;

Statistics getMeanMedianStd(
        const std::vector<  double > &i_vector
        )
{
   Statistics stats;
   std::vector <  double  > vector(i_vector);
   std::sort(vector.begin(),vector.end());
   double sumD = std::accumulate(vector.begin(),vector.end(),0.0);
   stats.mean = sumD/double(vector.size());
   std::vector < double  > Diff(vector.size());
   std::transform(vector.begin(),vector.end(),Diff.begin(),
                  std::bind2nd(std::minus< double  > (),stats.mean));
   stats.stdev  = std::inner_product(Diff.begin(),Diff.end(),
                                      Diff.begin(),0.0);
   stats.median = vector[std::floor(double(vector.size())/2.0)];
   stats.stdev = std::sqrt(stats.stdev/vector.size());
   return stats;
}

int main(int argc, char *argv[])
{
    std::vector < double  > values({10.0,12.0,11.0,4.5,10,13,22,2,1});
       for(unsigned int i=0; i < values.size(); ++i)
       {
          std::cout << values[i] << " " ;
       }
       std::cout << "\n";
       Statistics stats = getMeanMedianStd(values);
       std::cout << "Mean               = " << stats.mean   << "\n"
                 << "Median             = " << stats.median << "\n"
                 << "Standard Deviation = " << stats.stdev  << "\n";
    std::cout << "   ***   EXIT   ***\n";
    return EXIT_SUCCESS;
}

PS: for the median if the array contains an even number of values, then the bigger one is chosen between the two of them.

Tuesday, 26 July 2016

How to add metrics to DASOS

DASOS is our open source software for managing full-waveform (FW) LiDAR data [1] (http://miltomiltiadou.blogspot.co.uk/2015/03/las13vis.html).  Nevertheless it has a limited number of metrics, but it's design make it easy for people to add their own. So this blogpost aims to explain how to add your own metrics derived from the voxelised FW LiDAR data. It is advised to forward me (mmiltoo(at)gmail(dot)com) the new metrics classes in order to add the to future releases of DASOS and your name will also be added to the contributors.

There are two steps for generating your own metrics:
Step 1. Create a new Class for your new metric. To avoid confusion let's name our new class Metric. Please save the Metric.h file into the director ./include/Maps and the Metric.cpp file into the ./src/Maps directory to keep the files tidy. The Metrics class will inherit from the the base class Map. The header file (Metric.h) is standard and it should always be as the following (please replace all METRIC/Metric with the name of your actual new metric):

#ifndef METRIC_H
#define METRIC_H
#include "Map.h"
//-------------------------------------------------------------------------
/// @file Metric.h
/// @author <your name>
/// @version 1.0
/// @date <date generated>
/// @class Metric
/// @brief
//-------------------------------------------------------------------------


class Metric: public Map
{
public:
    //-------------------------------------------------------------------------
    /// @brief default constructor
    //-------------------------------------------------------------------------
    Metric(
            const std::string i_name,
            Volume *i_obj
            );
    //-------------------------------------------------------------------------
    /// @brief default destructor
    //-------------------------------------------------------------------------
    ~Metric();

private:
    //-------------------------------------------------------------------------
    /// @brief method that creates the Map
    //-------------------------------------------------------------------------
    void createMap();
};

#endif // METRIC_H


The only method that needs to be implemented is the createMap() which is a compulsory virtual function. The .cpp file should be as follow:

#include "Metric.h"

//-----------------------------------------------------------------------------
Metric::Metric(
        const std::string i_name,
        Volume *i_obj
        ):
    Map(i_name,i_obj)
{
}


//-----------------------------------------------------------------------------
void Metric::createMap()
{
   // Loop through all the voxels and generate the metric of interest
   // the variable m_noOfPixelsX, m_noOfPixelsY and m_noOfPixelsZ gives you 
   // the the number of voxels in x,y,z axis respectively
   for(unsigned int x=0; x<m_noOfPixelsX; ++x)
   {
      for(unsigned int y=0; y<m_noOfPixelsY; ++y)
      {
         for(unsigned int z=0; z<m_noOfPixelsZ; ++z)
         {
            // in m_mapValues all the values of the metrics are saved
            // the following command assigns the value 0 at the (x,y) position 
            // of the map
            m_mapValues[getIndex(x,y)]=-0.0f;

            // the voxelised 3D volume is the m_object variable and some 
            // useful examples of using it are shown below:

            // get the length of the voxel
            float voxelLength = m_object->getVoxelLen(); 
            // get the intesity at voxel (x,y,z)
            float intensity = m_object->getIntensity(x,y,z);
            // check whether the value of the voxel at (x,y,z) is considered to 
            // be inside or outside the scanned object. This checks whether the 
            // intensity value is above the boundary threshold. 
            bool isInside = m_object->isInside(x,y,z);            
         }
      }
   }
}

//-----------------------------------------------------------------------------
Metric::~Metric()
{}



Step 2: Link the new Metric Class with the rest of the program. This is done by modifying the .cpp file of MapsManager class.  Here it is shown the 4 additions that needs to be done in order to link your new metric with the rest of the program.


#include "MapsManager.h"
#include "ThicknessMap.h"
#include "FirstPatch.h"
#include "LastPatch.h"
// 1st ADDITION: include the header file of the new class
#include "Metric.h
// end of 1st ADDITION
#include <map>
#include <algorithm>

//-----------------------------------------------------------------------------
MapsManager::MapsManager():m_map(0),
    m_FWMetrics({"THICKNESS",
                 "LOWEST_RETURN",
                 "LAST_PATCH",
                 "FIRST_PATCH"
   // 2nd ADDITION: add the a name for your new metric
   // please do not forget the comma at the beginning
                 , "METRIC"
   // end of 2nd ADDITION
                })
{
   // The types should aggree with the fw metrics list
   m_types =
   {
      {"THICKNESS",3},
      {"FIRST_PATCH",9},
      {"LAST_PATCH",11}
    // 3rd ADDITION: give a number to your metric. This number must be unique
    // please do not forget the comma at the beginning
      , {"METRIC", 12}
   // end of 3rd ADDITION
   };
}

//-----------------------------------------------------------------------------
const std::vector<std::string> MapsManager::getNamesOfFWMetrics()const
{
   return m_FWMetrics;
}

//-----------------------------------------------------------------------------
void MapsManager::createMap(
        mapInfo *m_infoOfMap
        )
{
   if (m_map!=0)
   {
      delete m_map;
      m_map=0;
   }

   std::string s(m_infoOfMap->type);
   std::transform(s.begin(), s.end(), s.begin(), toupper);
   switch (m_types[s])
   {
  
   case 3:
      std::cout << "Density map\n";
      m_map = new DensityMap(m_infoOfMap->name,m_infoOfMap->obj);
      break;
  
  
   case 9:
      std::cout << "Length of first continues patch of non empty voxels\n";
      m_map = new FirstPatch(m_infoOfMap->name,m_infoOfMap->obj);
      break;
 
   case 11:
      std::cout << "Length of last continues patch of non empty voxels\n";
      m_map = new LastPatch(m_infoOfMap->name,m_infoOfMap->obj);
      break;

   // 4th ADDITION: link your metric class with the rest of the program
   case 12:
       std::cout << "Brief Discreption of the new Metric\n";
       m_map = new Metric(m_infoOfMap->name,m_infoOfMap->obj);
       break;
   // end of 4th ADDITION
   default:
      std::cout << std::string (s) << " is not a valid type of map";
      break;
   }
   // create and save the map
   if(m_map!=0)
   {
      m_map->createAndSave(m_infoOfMap->thres,m_infoOfMap->samp);
      delete m_map;
      m_map=0;
   }
}


//-----------------------------------------------------------------------------
MapsManager::~MapsManager()
{
   if(m_map!=0)
   {
      delete m_map;
   }
}


Once those steps are done you should be able to call your new metrics from the main program:

 ./DASOS -las myLasFile.LAS -map METRIC metric.asc


If you add the above class you should get a black asc file for all LAS files because all the values of the map are set to zero.

I hope you find that useful and if you have any questions please contact us at the following Google group:
https://groups.google.com/forum/#!forum/dasos---the-native-full-waveform-fw-lidar-software



Work Cited:
[1] Miltiadou, M., Grant, M. G., Campbell, N. D., Warren, M., Clewley, D., & Hadjimitsis, D. G. (2019, June). Open source software DASOS: Efficient accumulation, analysis, and visualisation of full-waveform lidar. In Seventh International Conference on Remote Sensing and Geoinformation of the Environment (RSCy2019) (Vol. 11174, p. 111741M). International Society for Optics and Photonics.

Friday, 29 April 2016

Bingo Game in C++

Once I was responsible of organising a bingo night and I wrote the following C++ code to make the game possible and allow myself to participate as well!

Since we were all programmers there, we enjoyed the idea of having our C++ code for playing the Bingo. By the end, I got slightly accused of customising the game in the favour of the event's organisers. Well, random is random and it wasn't my fault that we won the first two prizes!

The idea of the code is simple:
1. Firstly, an array of 100 elements is created and initialised to contained the values from 0 to 99 inclusive
2. Secondly, the array is shuttle using a random seed.
3. Then, the game starts! Every time the user presses 'enter' the next number inside the array is printed on the screen
4. and if there is a winner, the user may press 'b' to break the loop and end the game

The code of the Bingo is here:


#include < vector >
#include < iostream >
#include < time.h >
#include < stdlib.h >

int main(void)
{
    // create and initialise an array containing the values 0-99
    unsigned int legth = 100;
    std::vector < unsigned short int > m_bingoNumbers(legth);
    for(unsigned int i=0; i < m_bingoNumbers.size(); ++i)
    {
       m_bingoNumbers[i] = i;
    }

    // shuffle the array
    srand(time(NULL));
    for(unsigned int i=0; i< legth*2; ++i)
    {
       unsigned int rand1 = rand()%100;
       unsigned int rand2 = rand()%100;
       unsigned short int temp = m_bingoNumbers[rand1];
       m_bingoNumbers[rand1] = m_bingoNumbers[rand2];
       m_bingoNumbers[rand2] = temp;
    }

    // print numbers one by one 
    char b;
    for(unsigned int i=0; i< m_bingoNumbers.size(); ++i)
    {
       std::cin >> std::noskipws >> b;
       if(b=='b')
       {
           break;
       }
       if(b=='\n')
       {
          std::cout << " " << m_bingoNumbers[i];
       }
    }
    std::cout << "\n";
    return 0;
}


Here, there is also an example on how to compile and run the game:

$: g++ main.cpp -o bingo
$: ./bingo

 99
 4
 2
 50
 80
 41b



Monday, 29 February 2016

Finding all files of a given extension from a directory in C++

I have just found this useful recently and my blog it's a nice place to back it up for future reference. :)
The following small script takes as input an extension and a directory and prints all the files inside that directory that have that extension.


#include < iostream >
#include< stdio.h >
#include< cstdlib >
#include< iostream >
#include< string .h >
#include< fstream >
#include< sstream >
#include< dirent .h >
#include < vector >
// < extension > < directory >
int main(int argc, char *argv[])
{
    if(argc!=3)
    {
       std::cerr << "Too few arguments. Please include the following\n"
                 << "<.extension> \n";
       return EXIT_FAILURE;
    }
    std::string  dirStr(argv[2]);
    std::string extension(argv[1]);
    DIR *dir;
    struct dirent *ent;
    unsigned int count(0);
    if ((dir = opendir (dirStr.c_str())) != NULL) {
      /* print all the files and directories within directory */
      while ((ent = readdir (dir)) != NULL)
      {
        std::string current(ent->d_name);
        if(current.size()>extension.size())
        {
           std::string ext = current.substr(current.length()-extension.length());
           if(ext==extension)
           {
              std::cout << current << "\n";
              count++;
           }
        }
      }
    }
    std::cout << count << " files with extension " << extension << " found\n";
    std::cout << "   ***   EXIT   ***\n";
    return EXIT_SUCCESS;
}