ENH: add FSC_TOC and FSC_OG (top of canopy and on-ground fractional snow cover) (7a746782) · Commits · Remote Modules / let-it-snow

python/s2snow/snow_detector.py

+69 −16

Original line number	Diff line number	Diff line
		@@ -66,7 +66,7 @@ class snow_detector:
		# (if -1 the target resolution is equal to the max resolution of the input band)
		self.target_resolution = general.get("target_resolution", -1)

		# Parse vector option
		# Parse vector options
		vector_options = data["vector"]
		self.generate_vector = vector_options.get("generate_vector", True)
		self.generate_intermediate_vectors = vector_options.get("generate_intermediate_vectors", False)
		@@ -74,6 +74,15 @@ class snow_detector:
		self.gdal_trace_outline_dp_toler = vector_options.get("gdal_trace_outline_dp_toler", 0)
		self.gdal_trace_outline_min_area = vector_options.get("gdal_trace_outline_min_area", 0)

		# Parse FSC options
		if data.get("fsc", False):
		self.dofsc = data["fsc"]['dofsc']
		self.tcd_path = str(data["fsc"]['tcd'])
		self.fscOg_Eq = data["fsc"]['fscOg_Eq']
		self.fscToc_Eq = data["fsc"]['fscToc_Eq']
		else:
		self.dofsc = False

		# Parse cloud data
		cloud = data["cloud"]
		self.rf = cloud.get("rf")
		@@ -248,6 +257,9 @@ class snow_detector:
		self.composition_path = op.join(self.product_path, "LIS_COMPO.TIF")
		self.histogram_path = op.join(self.product_path, "LIS_HISTO.TXT")
		self.metadata_path = op.join(self.product_path, "LIS_METADATA.XML")
		self.ndsi_path = op.join(self.product_path, "LIS_NDSI.TIF")
		self.fscToc_path = op.join(self.product_path, "LIS_FSCTOC.TIF")
		self.fscOg_path = op.join(self.product_path, "LIS_FSCOG.TIF")

		def detect_snow(self, nbPass):
		# Set maximum ITK threads
		@@ -279,6 +291,8 @@ class snow_detector:
		self.pass1()
		if nbPass == 2:
		self.pass2()
		if self.dofsc:
		self.passfsc()

		# RGB composition
		composition_RGB(
		@@ -534,12 +548,12 @@ class snow_detector:
		logging.info("Start pass 1")

		# Pass1 : NDSI threshold
		ndsi_formula = "(im1b" + str(self.nGreen) + "-im1b" + str(self.nSWIR) + \
		self.ndsi_formula = "(im1b" + str(self.nGreen) + "-im1b" + str(self.nSWIR) + \
		")/(im1b" + str(self.nGreen) + "+im1b" + str(self.nSWIR) + ")"
		logging.info("ndsi formula: "+ ndsi_formula)
		logging.info("ndsi formula: "+ self.ndsi_formula)

		# NDSI condition (ndsi > x and not cloud)
		condition_ndsi = "(im2b1!=1 and (" + ndsi_formula + ")>" + str(self.ndsi_pass1) + " "
		condition_ndsi = "(im2b1!=1 and (" + self.ndsi_formula + ")>" + str(self.ndsi_pass1) + " "

		condition_pass1 = condition_ndsi + \
		" and im1b" + str(self.nRed) + "> " + str(self.rRed_pass1) + ")"
		@@ -708,12 +722,12 @@ class snow_detector:
		if snow_fraction > self.fsnow_total_lim:
		# Test zs value (-1 means that no zs elevation was found)
		if self.zs != -1:
		# NDSI threshold again
		ndsi_formula = "(im1b" + str(self.nGreen) + "-im1b" + str(self.nSWIR) + \
		")/(im1b" + str(self.nGreen) + "+im1b" + str(self.nSWIR) + ")"
		#~ # NDSI threshold again
		#~ ndsi_formula = "(im1b" + str(self.nGreen) + "-im1b" + str(self.nSWIR) + \
		#~ ")/(im1b" + str(self.nGreen) + "+im1b" + str(self.nSWIR) + ")"

		condition_pass2 = "(im3b1 != 1) and (im2b1>" + str(self.zs) + ")" \
		+ " and (" + ndsi_formula + "> " + str(self.ndsi_pass2) + ")" \
		+ " and (" + self.ndsi_formula + "> " + str(self.ndsi_pass2) + ")" \
		+ " and (im1b" + str(self.nRed) + ">" + str(self.rRed_pass2) + ")"

		bandMathPass2 = band_math([self.img,
		@@ -832,3 +846,42 @@ class snow_detector:
		self.ram,
		otb.ImagePixelType_uint8)
		bandMathPass3.ExecuteAndWriteOutput()
		bandMathPass3 = None

		def passfsc(self):
		# write NDSIx100 (0-100), nosnow (0) cloud (205) and nodata (254)
		expression = "(im2b1 == 100)?100*"+str(self.ndsi_formula)+":im2b1"
		bandMathApp = band_math([self.img,self.final_mask_path],
		self.ndsi_path,
		expression,
		self.ram,
		otb.ImagePixelType_uint8)
		bandMathApp.ExecuteAndWriteOutput()
		bandMathApp = None

		# write top-of-canopy FSC (0-100), nosnow (0) cloud (205) and nodata (254)
		#~ self.fscToc_Eq="1.45*ndsi-0.01"
		eq="min("+str(self.fscToc_Eq)+",1)"
		exp=eq.replace("ndsi", "im1b1/100") # ndsi was written in %
		expression = "(im2b1 == 100) ? 100*"+exp+" : im2b1"
		bandMathApp = band_math([self.ndsi_path,self.final_mask_path],
		self.fscToc_path,
		expression,
		self.ram,
		otb.ImagePixelType_uint8)
		bandMathApp.ExecuteAndWriteOutput()
		bandMathApp = None

		# write on-ground FSC (0-100), nosnow (0) cloud (205) and nodata (254)
		#~ self.fscOg_Eq="fscToc/(1-tcd)"
		eq="min("+str(self.fscOg_Eq)+",1)"
		exp=eq.replace("fscToc", "im1b1/100") # fscToc was written in %
		exp=exp.replace("tcd", "im3b1/100") # tcd is given in %
		expression = "(im2b1 == 100) ? 100*"+exp+" : im2b1"
		bandMathApp = band_math([self.fscToc_path,self.final_mask_path,self.tcd_path],
		self.fscOg_path,
		expression,
		self.ram,
		otb.ImagePixelType_uint8)
		bandMathApp.ExecuteAndWriteOutput()
		bandMathApp = None

styles/LIS_FSC_style_raster.qml

0 → 100644

+24 −0

Original line number	Diff line number	Diff line
		<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
		<qgis version="2.18.19" minimumScale="inf" maximumScale="1e+08" hasScaleBasedVisibilityFlag="0">
		<pipe>
		<rasterrenderer opacity="1" alphaBand="0" classificationMax="100" classificationMinMaxOrigin="User" band="1" classificationMin="0" type="singlebandpseudocolor">
		<rasterTransparency/>
		<rastershader>
		<colorrampshader colorRampType="DISCRETE" clip="0">
		<item alpha="255" value="0" label="No snow" color="#d6d2d0"/>
		<item alpha="255" value="20" label="0 - 20" color="#f0f9e8"/>
		<item alpha="255" value="40" label="20 - 40" color="#bae4bc"/>
		<item alpha="255" value="60" label="40 - 60" color="#7bccc4"/>
		<item alpha="255" value="80" label="60 - 80" color="#43a2ca"/>
		<item alpha="255" value="100" label="80 - 100" color="#0868ac"/>
		<item alpha="255" value="205" label="Cloud" color="#ffffff"/>
		<item alpha="255" value="inf" label="No data" color="#000000"/>
		</colorrampshader>
		</rastershader>
		</rasterrenderer>
		<brightnesscontrast brightness="0" contrast="0"/>
		<huesaturation colorizeGreen="128" colorizeOn="0" colorizeRed="255" colorizeBlue="128" grayscaleMode="0" saturation="0" colorizeStrength="100"/>
		<rasterresampler maxOversampling="2" zoomedOutResampler="bilinear"/>
		</pipe>
		<blendMode>0</blendMode>
		</qgis>

styles/LIS_NDSI_style_raster.qml

0 → 100644

+24 −0

Original line number	Diff line number	Diff line
		<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
		<qgis version="2.18.19" minimumScale="inf" maximumScale="1e+08" hasScaleBasedVisibilityFlag="0">
		<pipe>
		<rasterrenderer opacity="1" alphaBand="-1" classificationMax="100" classificationMinMaxOrigin="User" band="1" classificationMin="0" type="singlebandpseudocolor">
		<rasterTransparency/>
		<rastershader>
		<colorrampshader colorRampType="DISCRETE" clip="0">
		<item alpha="255" value="0" label="No snow" color="#d6d2d0"/>
		<item alpha="255" value="20" label="0 - 20" color="#f0f9e8"/>
		<item alpha="255" value="40" label="20 - 40" color="#bae4bc"/>
		<item alpha="255" value="60" label="40 - 60" color="#7bccc4"/>
		<item alpha="255" value="80" label="60 - 80" color="#43a2ca"/>
		<item alpha="255" value="100" label="80 - 100" color="#0868ac"/>
		<item alpha="255" value="205" label="Cloud" color="#ffffff"/>
		<item alpha="255" value="inf" label="No data" color="#000000"/>
		</colorrampshader>
		</rastershader>
		</rasterrenderer>
		<brightnesscontrast brightness="0" contrast="0"/>
		<huesaturation colorizeGreen="128" colorizeOn="0" colorizeRed="255" colorizeBlue="128" grayscaleMode="0" saturation="0" colorizeStrength="100"/>
		<rasterresampler maxOversampling="2" zoomedOutResampler="bilinear"/>
		</pipe>
		<blendMode>0</blendMode>
		</qgis>