Commit f1c3698d authored by Rémi Jugier's avatar Rémi Jugier
Browse files

water_mask options and cosims_mode

parent 1fd3d6cc
......@@ -55,8 +55,14 @@ conf_template = {"general":{"pout":"",
"dofsc": False,
"fscToc_Eq": "1.45*ndsi-0.01",
"fscOg_Eq": "fscToc/(1-tcd)",
"tcd": ""
}
"tcd": "",
"cosims_mode": False
},
"water_mask": {
"apply": False,
"path": None,
"raster_values": [1]
}
}
......@@ -304,7 +310,14 @@ def main():
group_cloud.add_argument("-strict_cloud_mask", type=str2bool, help="true/false")
group_fsc = parser.add_argument_group('fsc', 'fractional snow cover parameters')
group_fsc.add_argument("-fsc", type=str, help="path to tree cover density file, automatically activates sets fsc: dofsc to true", default='None')
group_fsc.add_argument("-fsc", type=str, help="path to tree cover density file, automatically activates sets fsc: dofsc to true")
group_fsc.add_argument("-cosims_mode", action='store_true', help="CoSIMS mode : Generate CoSIMS formatted outputs.")
group_water_mask = parser.add_argument_group('water_mask', 'water mask parameters')
group_water_mask.add_argument("-water_mask_path", type=str, help="Path to a raster or a shapefile")
group_water_mask.add_argument("-water_mask_raster_value", type=int, action='append', help="If the input water_mask_path is a raster, you can specify all the values corresponding " + \
"to water which are to be masked by repeating this optional argument -water_mask_raster_value value1 -water_mask_raster_value value2 etc... " + \
"If no values are specified, 1 will be used by default.")
args = parser.parse_args()
......@@ -416,11 +429,25 @@ def main():
logging.error("No DEM found!")
return 1
if args.fsc != 'None':
if args.fsc:
jsonData["fsc"]["dofsc"] = True
jsonData["fsc"]["tcd"] = args.fsc
jsonData["fsc"]["tcd"] = os.path.abspath(args.fsc)
jsonData["fsc"]["cosims_mode"] = args.cosims_mode
else:
jsonData["fsc"]["dofsc"] = False
if args.water_mask_path is not None:
jsonData["water_mask"]["apply"] = True
suffix = args.water_mask_path.split('.')[-1].lower()
if suffix not in ['shp', 'tif']:
raise IOError('Input water_mask_path must either be a GeoTIFF raster (.tif) or a shapefile (.shp)')
jsonData["water_mask"]["water_mask_path"] = os.path.abspath(args.water_mask_path)
if args.water_mask_raster_value is None:
jsonData["water_mask"]["water_mask_raster_values"] = [1]
else:
jsonData["water_mask"]["water_mask_raster_values"] = args.water_mask_raster_value
else:
jsonData["water_mask"]["apply"] = False
jsonFile = open(os.path.join(outputPath, "param_test.json"), "w")
jsonFile.write(json.dumps(jsonData, indent=4))
......
......@@ -36,7 +36,7 @@ from s2snow.app_wrappers import compute_snow_mask, compute_cloud_mask
from s2snow.app_wrappers import band_math, compute_snow_line
# Import utilities for snow detection
from s2snow.utils import polygonize, extract_band, burn_polygons_edges, composition_RGB
from s2snow.utils import polygonize, extract_band, burn_polygons_edges, composition_RGB, edit_raster_from_shapefile, edit_raster_from_raster, edit_nodata_value
from s2snow.utils import compute_percent, format_SEB_VEC_values, get_raster_as_array
# this allows GDAL to throw Python Exceptions
......@@ -80,9 +80,11 @@ class snow_detector:
self.tcd_path = str(data["fsc"]['tcd'])
self.fscOg_Eq = data["fsc"]['fscOg_Eq']
self.fscToc_Eq = data["fsc"]['fscToc_Eq']
self.cosims_mode = data["fsc"]['cosims_mode']
else:
self.dofsc = False
# Parse cloud data
cloud = data["cloud"]
self.rf = cloud.get("rf")
......@@ -133,6 +135,7 @@ class snow_detector:
gb_path_extracted = extract_band(inputs, "green_band", self.path_tmp, self.nodata)
rb_path_extracted = extract_band(inputs, "red_band", self.path_tmp, self.nodata)
sb_path_extracted = extract_band(inputs, "swir_band", self.path_tmp, self.nodata)
# Keep the input product directory basename as product_id
self.product_id = op.basename(op.dirname(inputs["green_band"]["path"]))
......@@ -196,6 +199,23 @@ class snow_detector:
yRes=self.target_resolution)
else:
sb_path_resampled = sb_path_extracted
#apply water mask as NAN values to extracted bands so that they are not used
if 'water_mask' in data:
if data['water_mask']['apply']:
water_mask_path = data['water_mask']['water_mask_path']
water_mask_type = water_mask_path.split('.')[-1].lower()
if water_mask_type == 'tif':
self.water_mask_raster_values = data['water_mask']['water_mask_raster_values']
edit_raster_from_raster(rb_path_resampled, water_mask_path, src_values=data['water_mask']['water_mask_raster_values'], applied_value=self.nodata)
edit_raster_from_raster(gb_path_resampled, water_mask_path, src_values=data['water_mask']['water_mask_raster_values'], applied_value=self.nodata)
edit_raster_from_raster(sb_path_resampled, water_mask_path, src_values=data['water_mask']['water_mask_raster_values'], applied_value=self.nodata)
elif water_mask_type == 'shp':
edit_raster_from_shapefile(rb_path_resampled, water_mask_path, applied_value=self.nodata)
edit_raster_from_shapefile(gb_path_resampled, water_mask_path, applied_value=self.nodata)
edit_raster_from_shapefile(sb_path_resampled, water_mask_path, applied_value=self.nodata)
else:
raise IOError('Input water_mask_path must either be a GeoTIFF raster (.tif) or a shapefile (.shp)')
# build vrt
logging.info("building bands vrt")
......@@ -232,7 +252,10 @@ class snow_detector:
self.label_no_snow = "0"
self.label_snow = "100"
self.label_cloud = "205"
self.label_no_data = "254"
if self.cosims_mode:
self.label_no_data = "255"
else:
self.label_no_data = "254"
# Build useful paths
self.pass1_path = op.join(self.path_tmp, "pass1.tif")
......@@ -294,42 +317,46 @@ class snow_detector:
if self.dofsc:
self.passfsc()
# RGB composition
composition_RGB(
self.img,
self.composition_path,
self.nSWIR,
self.nRed,
self.nGreen,
self.multi)
# Gdal polygonize (needed to produce composition)
# TODO: Study possible loss and issue with vectorization product
if self.generate_vector:
polygonize(
self.final_mask_path,
if self.cosims_mode:
self.create_cosims_metadata()
else:
# RGB composition
composition_RGB(
self.img,
self.composition_path,
self.nSWIR,
self.nRed,
self.nGreen,
self.multi)
# Gdal polygonize (needed to produce composition)
# TODO: Study possible loss and issue with vectorization product
if self.generate_vector:
polygonize(
self.final_mask_path,
self.final_mask_path,
self.final_mask_vec_path,
self.use_gdal_trace_outline,
self.gdal_trace_outline_min_area,
self.gdal_trace_outline_dp_toler)
# Burn polygons edges on the composition
# TODO add pass1 snow polygon in yellow
burn_polygons_edges(
self.composition_path,
self.final_mask_path,
self.final_mask_vec_path,
self.use_gdal_trace_outline,
self.gdal_trace_outline_min_area,
self.gdal_trace_outline_dp_toler)
# Burn polygons edges on the composition
# TODO add pass1 snow polygon in yellow
burn_polygons_edges(
self.composition_path,
self.final_mask_path,
self.label_snow,
self.label_cloud,
self.ram)
# Product formating
#~ format_SEB_VEC_values(self.final_mask_vec_path,
#~ self.label_snow,
#~ self.label_cloud,
#~ self.label_no_data)
self.create_metadata()
self.label_snow,
self.label_cloud,
self.ram)
# Product formating
#~ format_SEB_VEC_values(self.final_mask_vec_path,
#~ self.label_snow,
#~ self.label_cloud,
#~ self.label_no_data)
self.create_metadata()
def create_metadata(self):
# Compute and create the content for the product metadata file.
snow_percent = compute_percent(self.final_mask_path,
......@@ -356,6 +383,11 @@ class snow_detector:
name='CloudPercent').text = str(cloud_percent)
et = etree.ElementTree(root)
et.write(self.metadata_path, pretty_print=True)
def create_cosims_metadata(self):
self.create_metadata()
def extract_all_clouds(self):
if self.mode == 'lasrc':
......@@ -858,6 +890,8 @@ class snow_detector:
otb.ImagePixelType_uint8)
bandMathApp.ExecuteAndWriteOutput()
bandMathApp = None
if cosims_mode:
edit_nodata_value(self.ndsi_path, nodata_value=self.label_no_data)
# write top-of-canopy FSC (0-100), nosnow (0) cloud (205) and nodata (254)
#~ self.fscToc_Eq="1.45*ndsi-0.01"
......@@ -871,6 +905,7 @@ class snow_detector:
otb.ImagePixelType_uint8)
bandMathApp.ExecuteAndWriteOutput()
bandMathApp = None
edit_nodata_value(self.fscToc_path, nodata_value=self.label_no_data)
# write on-ground FSC (0-100), nosnow (0) cloud (205) and nodata (254)
#~ self.fscOg_Eq="fscToc/(1-tcd)"
......@@ -885,3 +920,4 @@ class snow_detector:
otb.ImagePixelType_uint8)
bandMathApp.ExecuteAndWriteOutput()
bandMathApp = None
edit_nodata_value(self.fscOg_path, nodata_value=self.label_no_data)
......@@ -14,6 +14,7 @@ from distutils import spawn
import numpy as np
import ogr
import gdal
import gdalconst
from gdalconst import GA_ReadOnly
......@@ -195,7 +196,6 @@ def extract_band(inputs, band, path_tmp, noData):
path = data_band["path"]
band_no = data_band["noBand"]
dataset = gdal.Open(path, GA_ReadOnly)
path_extracted = op.join(path_tmp, band+"_extracted.tif")
logging.info("extracting "+band)
......@@ -208,6 +208,60 @@ def extract_band(inputs, band, path_tmp, noData):
bandList=[band_no])
return path_extracted
def edit_nodata_value(raster_file, nodata_value=None, bands=None):
ds = gdal.Open(raster_file, gdal.GA_Update)
#iterate on each band
for band_no in range(1, ds.RasterCount+1):
if bands is not None:
if band_no not in bands:
#this band was not specified for edition, skip
continue
band = ds.GetRasterBand(band_no)
if nodata_value is None:
#remove nodata value
ds.GetRasterBand(band_no).DeleteNoDataValue()
else:
#change nodata value
ds.GetRasterBand(band_no).SetNoDataValue(nodata_value)
def edit_raster_from_shapefile(raster_target, src_shapefile, applied_value=0):
shape_mask = ogr.Open(src_shapefile)
ds = gdal.Open(raster_target, gdal.GA_Update)
for shape_mask_layer in shape_mask:
gdal.RasterizeLayer(ds, [1], shape_mask_layer, burn_values=[applied_value])
def edit_raster_from_raster(raster_target, src_raster, src_values, applied_value=0, layered_processing=False):
ds_mask = gdal.Open(src_raster, gdal.GA_ReadOnly)
band_mask = ds_mask.GetRasterBand(1)
ds = gdal.Open(raster_target, gdal.GA_Update)
band = ds.GetRasterBand(1)
if band.XSize != band_mask.XSize or band.YSize != band_mask.YSize:
raise IOError('array sizes from files do not match:\n%s'%('\n'.join([' - %s'%el for el in [raster_target, src_raster]])))
if layered_processing:
#iterate load line per line to avoid memory issues
for ii in range(band.YSize - 1, -1, -1):
data = band.ReadAsArray(xoff=0, yoff=ii, win_xsize=band.XSize, win_ysize=1, buf_xsize=band.XSize, buf_ysize=1)
data_mask = band_mask.ReadAsArray(xoff=0, yoff=ii, win_xsize=band.XSize, win_ysize=1, buf_xsize=band.XSize, buf_ysize=1)
for val in src_values:
data[data_mask==val] = applied_value
band.WriteArray(data, xoff=0, yoff=ii)
else:
data = band.ReadAsArray()
data_mask = band_mask.ReadAsArray()
for val in src_values:
data[data_mask==val] = applied_value
band.WriteArray(data)
def apply_color_table(raster_file_name, color_table):
......
Supports Markdown
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