Selection steep

execute.sh

@@ -5,8 +5,8 @@ INPUT_PATH=$2
 OUTPUT_PATH=$3
 METHOD=$4
 
-make clean
-make
+#make clean
+#make
 
 if [ $SENSOR == "LC08" ]
 then

scripts/collect-cpu-usage.sh

@@ -4,5 +4,5 @@ export LC_ALL="C"
 TIME_BETWEEN_COMMANDS=1
 echo TIMESTAMP, PID, USER, SYSTEM, GUEST, CPU, COMMAND
 while [ -e /proc/$1 ]; do
- /usr/bin/pidstat -p $1 -u $TIME_BETWEEN_COMMANDS 1 | grep -e '[0-9]\{2\}:[0-9]\{2\}:[0-9]\{2\}[ \t]*[0-9]' | awk -v date="$( date +"%s" )" '{print date", "$3", "$4", "$5", "$6", "$7", "$9; }' 2> /dev/null
+ /usr/bin/pidstat -p $1 -u -l $TIME_BETWEEN_COMMANDS 1 | grep -e '[0-9]\{2\}:[0-9]\{2\}:[0-9]\{2\}[ \t]*[0-9]' | awk -v date="$( date +"%s" )" '{print date", "$3", "$4", "$5", "$6", "$7", "$9; }' 2> /dev/null
 done

src/STEEP.cpp

@@ -0,0 +1,153 @@
+#include "STEEP.h"
+
+void compute_H0(double net_radiation_line[], double soil_heat_flux[], int width_band, double ho_line[])
+{
+ for (int col = 0; col < width_band; col++)
+ ho_line[col] = net_radiation_line[col] - soil_heat_flux[col];
+};
+
+void filter_valid_values(double *target_line, double *target_values, int width_band, int *pos)
+{
+ for (int col = 0; col < width_band; col++)
+ {
+ if (!isnan(target_line[col]) && !isinf(target_line[col]))
+ {
+ target_values[*pos] = target_line[col];
+ (*pos)++;
+ }
+ }
+}
+
+void get_quartiles(TIFF *target, double *v_quartile, int height_band, int width_band, double first_interval, double last_interval)
+{
+ const int SIZE = height_band * width_band;
+ double target_line[width_band];
+ double *target_values = (double *)malloc(sizeof(double) * SIZE);
+
+ if (target_values == NULL)
+ exit(15);
+
+ int pos = 0;
+ for (int line = 0; line < height_band; line++)
+ {
+ read_line_tiff(target, target_line, line);
+ filter_valid_values(target_line, target_values, width_band, &pos);
+ }
+
+ sort(target_values, target_values + pos);
+
+ v_quartile[0] = target_values[int(floor(first_interval * pos))];
+ v_quartile[1] = target_values[int(floor(last_interval * pos))];
+
+ free(target_values);
+}
+
+Candidate getHotPixelSTEPP(TIFF** ndvi, TIFF** surface_temperature, TIFF** albedo, TIFF** net_radiation, TIFF** soil_heat, int height_band, int width_band)
+{
+ std::vector<double> ndvi_line(width_band);
+ std::vector<double> surface_temperature_line(width_band);
+ std::vector<double> net_radiation_line(width_band);
+ std::vector<double> soil_heat_line(width_band);
+ std::vector<double> albedo_line(width_band);
+ std::vector<double> ho_line(width_band);
+ std::vector<Candidate> candidatesGroup;
+
+ std::vector<double> ndviQuartile(2);
+ std::vector<double> tsQuartile(2);
+ std::vector<double> albedoQuartile(2);
+
+ get_quartiles(*ndvi, ndviQuartile.data(), height_band, width_band, 0.15, 0.85);
+ get_quartiles(*albedo, albedoQuartile.data(), height_band, width_band, 0.50, 0.75);
+ get_quartiles(*surface_temperature, tsQuartile.data(), height_band, width_band, 0.85, 0.97);
+
+ for (int line = 0; line < height_band; line++)
+ {
+ read_line_tiff(*net_radiation, net_radiation_line.data(), line);
+ read_line_tiff(*soil_heat, soil_heat_line.data(), line);
+
+ compute_H0(net_radiation_line.data(), soil_heat_line.data(), width_band, ho_line.data());
+
+ read_line_tiff(*ndvi, ndvi_line.data(), line);
+ read_line_tiff(*surface_temperature, surface_temperature_line.data(), line);
+ read_line_tiff(*albedo, albedo_line.data(), line);
+
+ for (int col = 0; col < width_band; col++)
+ {
+ bool ndviValid = !std::isnan(ndvi_line[col]) && ndvi_line[col] > 0.10 && ndvi_line[col] < ndviQuartile[0];
+ bool albedoValid = !std::isnan(albedo_line[col]) && albedo_line[col] > albedoQuartile[0] && albedo_line[col] < albedoQuartile[1];
+ bool tsValid = !std::isnan(surface_temperature_line[col]) && surface_temperature_line[col] > tsQuartile[0] && surface_temperature_line[col] < tsQuartile[1];
+
+ if (albedoValid && ndviValid && tsValid)
+ {
+ candidatesGroup.emplace_back(ndvi_line[col], surface_temperature_line[col], net_radiation_line[col],
+ soil_heat_line[col], ho_line[col], line, col);
+ }
+ }
+ }
+
+ if (candidatesGroup.empty())
+ {
+ cerr << "Pixel problem! - There are no final candidates";
+ exit(15);
+ }
+
+ std::sort(candidatesGroup.begin(), candidatesGroup.end(), compare_candidate_temperature);
+ unsigned int pos = static_cast<unsigned int>(std::floor(candidatesGroup.size() * 0.5));
+
+ return candidatesGroup[pos];
+}
+
+Candidate getColdPixelSTEPP(TIFF** ndvi, TIFF** surface_temperature, TIFF** albedo, TIFF** net_radiation, TIFF** soil_heat, int height_band, int width_band)
+{
+ std::vector<double> ndvi_line(width_band);
+ std::vector<double> surface_temperature_line(width_band);
+ std::vector<double> net_radiation_line(width_band);
+ std::vector<double> soil_heat_line(width_band);
+ std::vector<double> albedo_line(width_band);
+ std::vector<double> ho_line(width_band);
+ std::vector<Candidate> candidatesGroup;
+
+ std::vector<double> ndviQuartile(2);
+ std::vector<double> tsQuartile(2);
+ std::vector<double> albedoQuartile(2);
+
+ get_quartiles(*ndvi, ndviQuartile.data(), height_band, width_band, 0.15, 0.97);
+ get_quartiles(*albedo, albedoQuartile.data(), height_band, width_band, 0.25, 0.50);
+ get_quartiles(*surface_temperature, tsQuartile.data(), height_band, width_band, 0.20, 0.85);
+
+ for (int line = 0; line < height_band; line++)
+ {
+ read_line_tiff(*net_radiation, net_radiation_line.data(), line);
+ read_line_tiff(*soil_heat, soil_heat_line.data(), line);
+
+ compute_H0(net_radiation_line.data(), soil_heat_line.data(), width_band, ho_line.data());
+
+ read_line_tiff(*ndvi, ndvi_line.data(), line);
+ read_line_tiff(*surface_temperature, surface_temperature_line.data(), line);
+ read_line_tiff(*albedo, albedo_line.data(), line);
+
+ for (int col = 0; col < width_band; col++)
+ {
+ bool ndviValid = !std::isnan(ndvi_line[col]) && ndvi_line[col] > ndviQuartile[1];
+ bool albedoValid = !std::isnan(albedo_line[col]) && albedo_line[col] > albedoQuartile[0] && albedo_line[col] < albedoQuartile[1];
+ bool tsValid = !std::isnan(surface_temperature_line[col]) && surface_temperature_line[col] < tsQuartile[0];
+
+ if (albedoValid && ndviValid && tsValid)
+ {
+ candidatesGroup.emplace_back(ndvi_line[col], surface_temperature_line[col], net_radiation_line[col],
+ soil_heat_line[col], ho_line[col], line, col);
+ }
+ }
+ }
+
+ if (candidatesGroup.empty())
+ {
+ cerr << "Pixel problem! - There are no final candidates";
+ exit(15);
+ }
+
+ std::sort(candidatesGroup.begin(), candidatesGroup.end(), compare_candidate_temperature);
+ unsigned int pos = static_cast<unsigned int>(std::floor(candidatesGroup.size() * 0.5));
+
+ return candidatesGroup[pos];
+}

src/STEEP.h

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "types.h"
+#include "utils.h"
+#include "parameters.h"
+#include "pixel_reader.h"
+#include "candidate.h"
+
+void compute_H0(double net_radiation_line[], double soil_heat_flux[], int width_band, double ho_line[]);
+
+void filter_valid_values(double *target_line, double *target_values, int width_band, int *pos);
+
+void get_quartiles(TIFF *target, double *v_quartile, int height_band, int width_band, double first_interval, double last_interval);
+
+Candidate getHotPixelSTEPP(TIFF** ndvi, TIFF** surface_temperature, TIFF** albedo, TIFF** net_radiation, TIFF** soil_heat, int height_band, int width_band);
+
+Candidate getColdPixelSTEPP(TIFF** ndvi, TIFF** surface_temperature, TIFF** albedo, TIFF** net_radiation, TIFF** soil_heat, int height_band, int width_band);

src/asebal.cpp

@@ -117,20 +117,30 @@ Candidate getHotPixel(TIFF** ndvi, TIFF** surface_temperature, TIFF** albedo, TI
  }
 
  }
+
+ if(candidatesGroupI.size() <= 0) {
+ cerr << "Pixel problem! - There are no precandidates";
+ exit(15);
+ }
 
  //Creating second pixel group, all values lower than the 3rd quartile are excluded
  sort(candidatesGroupI.begin(), candidatesGroupI.end(), compare_candidate_temperature);
  unsigned int pos = int(floor(candidatesGroupI.size() * 0.75));
  vector<Candidate> candidatesGroupII(candidatesGroupI.begin() + pos, candidatesGroupI.end());
 
+ if(candidatesGroupII.size() <= 0) {
+ cerr << "Pixel problem! - There are no final candidates";
+ exit(15);
+ }
+
  pos = int(floor(candidatesGroupII.size() * 0.5));
  Candidate hotPixel = candidatesGroupII[pos];
 
  free(ndviQuartile);
  free(tsQuartile);
  free(albedoQuartile);
 
- hotPixel.toString();
+ //hotPixel.toString();
 
  return hotPixel;
 }
@@ -193,19 +203,31 @@ Candidate getColdPixel(TIFF** ndvi, TIFF** surface_temperature, TIFF** albedo, T
 
  }
 
+ if(candidatesGroupI.size() <= 0) {
+ cerr << "Pixel problem! - There are no precandidates";
+ exit(15);
+ }
+
  //Creating second pixel group, all pixels with values higher than the 1 st quartile are excluded
  sort(candidatesGroupI.begin(), candidatesGroupI.end(), compare_candidate_temperature);
  unsigned int pos = int(floor(candidatesGroupI.size() * 0.25));
  vector<Candidate> candidatesGroupII(candidatesGroupI.begin(), candidatesGroupI.begin() + pos);
+
 
+ if(candidatesGroupII.size() <= 0) {
+ cerr << "Pixel problem! - There are no final candidates";
+ exit(15);
+ }
+
  pos = int(floor(candidatesGroupII.size() * 0.5));
+
  Candidate coldPixel = candidatesGroupII[pos];
 
  free(ndviQuartile);
  free(tsQuartile);
  free(albedoQuartile);
 
- coldPixel.toString();
+ // coldPixel.toString();
 
  return coldPixel;
-}
+}

src/esaSel.cpp

@@ -503,6 +503,11 @@ pair<Candidate, Candidate> esaPixelSelect(TIFF** ndvi, TIFF** surface_temperatur
 
  }
 
+ if(listTS.size() <= 0) {
+ cerr << "Pixel problem! - There are no precandidates";
+ exit(15);
+ }
+
  vector<Candidate> listNDVI (listTS);
 
  sort(listTS.begin(), listTS.end(), compare_candidate_temperature);
@@ -533,6 +538,11 @@ pair<Candidate, Candidate> esaPixelSelect(TIFF** ndvi, TIFF** surface_temperatur
 
  }
 
+ if(final_histTS.size() <= 0) {
+ cerr << "Pixel problem! - There are no final TS candidates";
+ exit(15);
+ }
+
  vector<Candidate> histNDVI[binNDVI], final_histNDVI;
  for(Candidate c : listNDVI) {
 
@@ -552,6 +562,11 @@ pair<Candidate, Candidate> esaPixelSelect(TIFF** ndvi, TIFF** surface_temperatur
 
  }
 
+ if(final_histNDVI.size() <= 0) {
+ cerr << "Pixel problem! - There are no final NDVI candidates";
+ exit(15);
+ }
+
  // Select cold pixel
  int pixel_count = 0, n1 = 1, n2 = 1, ts_pos, ndvi_pos, beginTs = 0, beginNDVI = final_histNDVI.size() - 1;
  vector<Candidate> coldPixels;
@@ -586,6 +601,11 @@ pair<Candidate, Candidate> esaPixelSelect(TIFF** ndvi, TIFF** surface_temperatur
 
  }
 
+ if(coldPixels.size() <= 0) {
+ cerr << "Pixel problem! - There are no cold candidates";
+ exit(15);
+ }
+
  //Select hot pixel
  pixel_count = 0, n1 = 1, n2 = 1;
  vector<Candidate> hotPixels;
@@ -621,6 +641,11 @@ pair<Candidate, Candidate> esaPixelSelect(TIFF** ndvi, TIFF** surface_temperatur
 
  }
 
+ if(hotPixels.size() <= 0) {
+ cerr << "Pixel problem! - There are no hot candidates";
+ exit(15);
+ }
+
  sort(coldPixels.begin(), coldPixels.end(), compare_candidate_ndvi);
  sort(hotPixels.begin(), hotPixels.end(), compare_candidate_temperature);
 

src/landsat.cpp

@@ -79,8 +79,11 @@ void Landsat::process_partial_products(TIFF* read_bands[], MTL mtl, Station stat
  double eo_emissivity_line[width_band], ea_emissivity_line[width_band], enb_emissivity_line[width_band];
  double large_wave_radiation_atmosphere_line[width_band], large_wave_radiation_surface_line[width_band];
  double short_wave_radiation_line[width_band];
-
+ 
  //Calculating the partial products for each line
+
+ //DEBUG - OLHAR ESSES CARAS AQUI!!
+
  for(int line = 0; line < height_band; line ++){
  radiance_function(read_bands, mtl, sensor, width_band, line, radiance_line, this->noData);
  reflectance_function(read_bands, mtl, sensor, radiance_line, width_band, line, reflectance_line, this->noData);
@@ -141,16 +144,15 @@ void Landsat::process_final_products(Station station, MTL mtl){
 
  Candidate hot_pixel, cold_pixel;
 
- if (this->method == 0) { //Our SEBAL
- hot_pixel = select_hot_pixel(&ndvi, &surface_temperature, &net_radiation, &soil_heat, height_band, width_band);
- cold_pixel = select_cold_pixel(&ndvi, &surface_temperature, &net_radiation, &soil_heat, height_band, width_band);
+ if (this->method == 0) { //STEEP
+ hot_pixel = getHotPixelSTEPP(&ndvi, &surface_temperature, &albedo, &net_radiation, &soil_heat, height_band, width_band);
+ cold_pixel = getColdPixelSTEPP(&ndvi, &surface_temperature, &albedo, &net_radiation, &soil_heat, height_band, width_band);
  } else if (this->method == 1) { //ASEBAL
  hot_pixel = getHotPixel(&ndvi, &surface_temperature, &albedo, &net_radiation, &soil_heat, height_band, width_band);
  cold_pixel = getColdPixel(&ndvi, &surface_temperature, &albedo, &net_radiation, &soil_heat, height_band, width_band);
  } else if (this->method == 2) { //ESA SEBAL
  TIFF* land_cover = TIFFOpen(this->land_cover_path.c_str(), "r");
  pair<Candidate, Candidate> pixels = esaPixelSelect(&ndvi, &surface_temperature, &albedo, &net_radiation, &soil_heat, &land_cover, height_band, width_band, this->output_path);
-
  hot_pixel = pixels.first, cold_pixel = pixels.second;
  }
 

src/landsat.h

@@ -6,6 +6,7 @@
 #include "parameters.h"
 #include "asebal.h"
 #include "esaSel.h"
+#include "STEEP.h"
 
 /**
  * @brief Struct to manage the products calculation.

src/main.cpp

@@ -54,7 +54,7 @@ int main(int argc, char *argv[]){
  if(noData_flag.substr(0,5) == "-nan=")
  noData = atof(noData_flag.substr(5, noData_flag.size()).c_str());
  }
-
+  
  string tal_path = argv[9];
 
  TIFF *bands_resampled[8];
@@ -70,11 +70,11 @@ int main(int argc, char *argv[]){
 
  Landsat landsat = Landsat(tal_path, output_path, method, noData, landCoverPath);
 
- begin = chrono::steady_clock::now();
- landsat.process_partial_products(bands_resampled, mtl, station, sensor);
- end = chrono::steady_clock::now();
- time_span_us = chrono::duration_cast< chrono::duration<double, micro> >(end - begin);
- printf("PHASE 1 - DURATION, %.5f\n", time_span_us);
+ // begin = chrono::steady_clock::now();
+ // landsat.process_partial_products(bands_resampled, mtl, station, sensor);
+ // end = chrono::steady_clock::now();
+ // time_span_us = chrono::duration_cast< chrono::duration<double, micro> >(end - begin);
+ // printf("PHASE 1 - DURATION, %.5f\n", time_span_us);
 
  begin = chrono::steady_clock::now();
 
@@ -84,4 +84,4 @@ int main(int argc, char *argv[]){
  time_span_us = chrono::duration_cast< chrono::duration<double, micro> >(end - begin);
  printf("PHASE 2 - DURATION, %.5f\n", time_span_us);
  return 0;
-}
+}