Commit 84840b4f authored by Manuel Grizonnet's avatar Manuel Grizonnet
Browse files

ENH: fix nosnow value in histogram;compute_snowline executable and snow_line...

ENH: fix nosnow value in histogram;compute_snowline executable and snow_line utils to investigate cold clouds
parent df12420a
#!/usr/bin/env python
import sys
import os.path as op
import json
import argparse
def show_help():
"""Show help of the build_json script for theia N2A products"""
print "This script is used to build json configuration file use then to compute snow mask using OTB applications on Spot/LandSat/Sentinel-2 products from theia platform"
print "Usage: python build_theia_json -s [landsat|s2|take5] -d image_directory -e srtm_tile -o file.json"
print "python run_snow_detector.py help to show help"
#----------------- MAIN ---------------------------------------------------
def main():
""" Script to build json from theia N2A product"""
parser = argparse.ArgumentParser(description='Build json from THEIA product')
parser.add_argument("-s", help="select input sensors")
parser.add_argument("-d", help="input dir")
parser.add_argument("-o", help="input dir")
parser.add_argument("-do", help="input dir")
args = parser.parse_args()
#print(args.accumulate(args.integers))
#Parse sensor
if (args.s == 's2'):
multi=10
#Build json file
data = {}
data["general"]={
"pout":args.do,
"nodata":-10000,
"ram":1024,
"nb_threads":1,
"generate_vector":"false",
"preprocessing":"false",
"log":"true",
"multi":10
}
data["cloud"]={
"shadow_mask":32,
"all_cloud_mask":1,
"high_cloud_mask":128,
"rf":12,
"red_darkcloud":500,
"red_backtocloud":100
}
data["snow"]={
"dz":100,
"ndsi_pass1":0.4,
"red_pass1":200,
"ndsi_pass2":0.15,
"red_pass2":120,
"fsnow_lim":0.1,
"fsnow_total_lim":0.001
}
fp = open(args.o, 'w')
fp.write(json.dumps(data,indent=4, sort_keys=True))
fp.close()
if __name__ == "__main__":
main()
......@@ -353,16 +353,17 @@ class snow_detector :
#Pass 2: compute snow fraction (c++)
nb_snow_pixels = histo_utils_ext.compute_nb_pixels_between_bounds(self.ndsi_pass1_path, 0 , 255)
print "Number of snow pixels ", nb_snow_pixels
if (nb_snow_pixels > self.fsnow_total_lim):
#Pass 2: determine the Zs elevation fraction (c++)
#Save histogram values for logging
histo_log=op.join(self.path_tmp,"histogram.txt")
#c++ function
self.zs=histo_utils_ext.compute_snowline(self.dem,self.ndsi_pass1_path,op.join(self.path_tmp,"cloud_pass1.tif"), self.dz, self.fsnow_lim, False, -2, -self.dz/2, histo_log)
print "computed ZS:", self.zs
#Compute Zs elevation fraction and histogram values
#We compute it in all case as we need to check histogram values to detect cold clouds in optionnal pass4
histo_log=op.join(self.path_tmp,"histogram.txt")
#c++ function
self.zs=histo_utils_ext.compute_snowline(self.dem,self.ndsi_pass1_path,op.join(self.path_tmp,"cloud_pass1.tif"), self.dz, self.fsnow_lim, False, -2, -self.dz/2, histo_log)
print "computed ZS:", self.zs
if (nb_snow_pixels > self.fsnow_total_lim):
#Test zs value (-1 means that no zs elevation was found)
if (self.zs !=-1):
#NDSI threshold again
......
......@@ -17,3 +17,7 @@ target_link_libraries(compute_snow_mask ${OTB_LIBRARIES})
#Build cloud mask exe
add_executable(compute_cloud_mask compute_cloud_mask.cxx)
target_link_libraries(compute_cloud_mask ${OTB_LIBRARIES})
#Build snow line exe
add_executable(compute_snowline compute_snowline.cxx)
target_link_libraries(compute_snowline histo_utils)
#include "histo_utils.h"
int main(int argc, char * argv[])
{
if (argc != 10)
{
std::cout << "infname inmasksnowfname inmaskcloudfname dz fsnow_lim reverse offset center_offset histo_file" << std::endl;
}
return compute_snowline(argv[1], argv[2], argv[3], atoi(argv[4]), atof(argv[5]), atoi(argv[6]), atoi(argv[7]), atoi(argv[8]), argv[9]);
}
......@@ -226,11 +226,12 @@ void print_histogram (const itk::Statistics::ImageToHistogramFilter<itk::VectorI
idx4[2] = 1;
const HistogramType::AbsoluteFrequencyType z_center = histogram.GetMeasurementVector(idx1)[0];
const int Nz = histogram.GetFrequency(idx1) + histogram.GetFrequency(idx2);
const int Nz = histogram.GetFrequency(idx1) + histogram.GetFrequency(idx2) + histogram.GetFrequency(idx3) + histogram.GetFrequency(idx4);
const int fcloud_z = histogram.GetFrequency(idx3) + histogram.GetFrequency(idx4);
const int fsnow_z = histogram.GetFrequency(idx2) + histogram.GetFrequency(idx4);
const int fnosnow_z = Nz - (fcloud_z + fsnow_z);
const int fnosnow_z = histogram.GetFrequency(idx1);
//myfile << histogram.GetFrequency(idx1) << ", " << histogram.GetFrequency(idx2) << ", "<< histogram.GetFrequency(idx3) << ", "<< histogram.GetFrequency(idx4) << std::endl;
myfile << z_center << ", " << Nz << ", "<< fcloud_z << ", "<< fsnow_z << ", "<< fnosnow_z << std::endl;
}
......
#!/usr/bin/python
#coding=utf8
#=========================================================================
#
# Program: lis
# Language: Python
#
# Copyright (c) Simon Gascoin
# Copyright (c) Manuel Grizonnet
#
# See lis-copyright.txt for details.
#
# This software is distributed WITHOUT ANY WARRANTY; without even
# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the above copyright notices for more information.
#
#=========================================================================
import numpy as np
from matplotlib import pyplot as plt
from optparse import OptionParser
def load_histo(histo_path):
with open(histo_path) as f:
v = np.loadtxt(f, delimiter=",", dtype='float', comments="#", skiprows=3, usecols=(0,3,4))
fsnow_rate=v[:,1]/(v[:,1]+v[:,2])
#print v[:,1]
#print v[:,2]
#b = np.zeros(6).reshape(2, 3)
#print fsnow_rate
print fsnow_rate[0]
print fsnow_rate
print np.shape(fsnow_rate)[0]
plt.plot(np.arange(np.shape(fsnow_rate)[0]), fsnow_rate[:], 'ro')
#plt.axis([0, 6, 0, 20])
plt.show()
def print_histo(histo_path):
with open(histo_path) as f:
v = np.loadtxt(f, delimiter=",", dtype='int', comments="#", skiprows=3, usecols=(0,1,3))
#_hist = np.ravel(v) # 'flatten' v
#fig = plt.figure()
#ax1 = fig.add_subplot(111)
#n, bins, patches = ax1.hist(v_hist, bins=50, normed=1, facecolor='green')
#plt.show()
print v
dem=v[:,0]
width = 0.8
#indices = np.arange(len(dem))
plt.bar(dem, v[:,1], width=width, color="red", label="all")
plt.bar([i+0.25*width for i in dem], v[:,2], width=0.5*width, color="blue", alpha=1. , label="snow")
plt.xticks(dem+width/2.,
[i for i in dem] )
plt.legend()
plt.show()
def main():
parser = OptionParser(usage="usage: %prog [options] filename",
version="%prog")
parser.add_option("-f","--file", help="absolute path to histogram file", dest="histo_path", type="string")
#parser.add_option("-f","--file", help="absolute path to histogram file", dest="histo_path", type="string")
(opts, args) = parser.parse_args()
if opts.histo_path is None:
print "A mandatory option is missing\n"
parser.print_help()
exit(-1)
else:
#print opts.path
#print_histo(opts.histo_path)
load_histo(opts.histo_path)
if __name__ == '__main__':
main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment