Commit 0e0fceea authored by Manuel Grizonnet's avatar Manuel Grizonnet
Browse files

ENH: run autopep8 on python scripts to improve code quality

parent de68760d
......@@ -5,7 +5,8 @@ import os.path as op
import json
from s2snow import snow_detector
VERSION="1.2"
VERSION = "1.2"
def show_help():
"""Show help of the run_snow_detector script"""
......@@ -14,34 +15,35 @@ def show_help():
print "python run_snow_detector.py version to show version"
print "python run_snow_detector.py help to show help"
def show_version():
print VERSION
#----------------- MAIN ---------------------------------------------------
# ----------------- MAIN ---------------------------------------------------
def main(argv):
""" main script of snow extraction procedure"""
json_file=argv[1]
json_file = argv[1]
#load json_file from json files
# Load json_file from json files
with open(json_file) as json_data_file:
data = json.load(json_data_file)
data = json.load(json_data_file)
general = data["general"]
pout = general.get("pout")
log = general.get("log", True)
if log:
sys.stdout = open(op.join(pout, "stdout.log"), 'w')
sys.stderr = open(op.join(pout, "stderr.log"), 'w')
sd = snow_detector.snow_detector(data)
sd.detect_snow(2)
if __name__ == "__main__":
if len(sys.argv) != 2 :
if len(sys.argv) != 2:
show_help()
else:
if sys.argv[1] == "version":
......
import sys, os, numpy, random
import sys
import os
import numpy
import random
import os.path as op
import gdal, gdalconst
import gdal
import gdalconst
from subprocess import call
def show_help():
print "This script is used to create clouds on data"
print "Usage: cloud_builder.py mode plaincloudthreshold randomcloudthreshold inputpath outputplaincloudpath ouputrandomcloudpath"
print "Mode : 0 %plain cloud image, 1 %random cloud image, 2 both"
print "This script is used to create clouds on data"
print "Usage: cloud_builder.py mode plaincloudthreshold randomcloudthreshold inputpath outputplaincloudpath ouputrandomcloudpath"
print "Mode : 0 %plain cloud image, 1 %random cloud image, 2 both"
def main(argv):
mode = argv[1]
mode = int(mode)
......@@ -17,35 +24,38 @@ def main(argv):
input_path = argv[4]
output_plain_cloud_path = argv[5]
output_random_cloud_path = argv[6]
dataset = gdal.Open(input_path, gdalconst.GA_ReadOnly)
wide = dataset.RasterXSize
high = dataset.RasterYSize
if(mode == 0 or mode == 2):
#build half cloud image
str_exp = "idxX>"+str(int(wide*plain_cloud_threshold))+"?im1b1:205"
call(["otbcli_BandMathX", "-il", input_path, "-out", output_plain_cloud_path, "-exp", str_exp])
# build half cloud image
str_exp = "idxX>" + \
str(int(wide * plain_cloud_threshold)) + "?im1b1:205"
call(["otbcli_BandMathX", "-il", input_path, "-out",
output_plain_cloud_path, "-exp", str_exp])
if(mode == 1 or mode == 2):
#build random cloud image
# build random cloud image
band = dataset.GetRasterBand(1)
array = band.ReadAsArray(0, 0, wide, high)
for i in range(0, wide):
for j in range(0, high):
if(random.randint(0, 100) < random_cloud_threshold):
array[i,j] = 205
output_random_cloud_raster = gdal.GetDriverByName('GTiff').Create(output_random_cloud_path, wide, high, 1 ,gdal.GDT_Byte)
output_random_cloud_raster.GetRasterBand(1).WriteArray(array)
array[i, j] = 205
output_random_cloud_raster = gdal.GetDriverByName('GTiff').Create(
output_random_cloud_path, wide, high, 1, gdal.GDT_Byte)
output_random_cloud_raster.GetRasterBand(1).WriteArray(array)
output_random_cloud_raster.FlushCache()
# georeference the image and set the projection
output_random_cloud_raster.SetGeoTransform(dataset.GetGeoTransform())
output_random_cloud_raster.SetProjection(dataset.GetProjection())
output_random_cloud_raster.SetProjection(dataset.GetProjection())
if __name__ == "__main__":
if len(sys.argv) != 7:
show_help()
else:
main(sys.argv)
This diff is collapsed.
......@@ -5,25 +5,30 @@ import sys
import subprocess
import ast
def show_help():
print "This script is used to compute srtm mask from a vrt file to a region extent"
print "Usage: preprocessing.py srtm.vrt img.tif output.tif"
print "This script is used to compute srtm mask from a vrt file to a region extent"
print "Usage: preprocessing.py srtm.vrt img.tif output.tif"
def get_extent(geotransform, cols, rows):
extent=[]
xarr=[0, cols]
yarr=[0, rows]
extent = []
xarr = [0, cols]
yarr = [0, rows]
for px in xarr:
for py in yarr:
x=geotransform[0]+(px*geotransform[1])+(py*geotransform[2])
y=geotransform[3]+(px*geotransform[4])+(py*geotransform[5])
x = geotransform[0] + (px * geotransform[1]) + \
(py * geotransform[2])
y = geotransform[3] + (px * geotransform[4]) + \
(py * geotransform[5])
extent.append([x, y])
yarr.reverse()
return extent
def build_dem(psrtm, pimg, pout):
#load datasets
# load datasets
source_dataset = gdal.Open(psrtm, gdalconst.GA_Update)
source_geotransform = source_dataset.GetGeoTransform()
source_projection = source_dataset.GetProjection()
......@@ -34,37 +39,52 @@ def build_dem(psrtm, pimg, pout):
wide = target_dataset.RasterXSize
high = target_dataset.RasterYSize
#compute extent xminymin and yminymax
# compute extent xminymin and yminymax
extent = get_extent(target_geotransform, wide, high)
print("Extent: " + str(extent))
te = "".join(str(extent[1] + extent[3])) #xminymin xmaxymax
te = "".join(str(extent[1] + extent[3])) # xminymin xmaxymax
te = ast.literal_eval(te)
te = ' '.join([str(x) for x in te])
print(te)
#get target resolution
resolution = target_geotransform[1] #or geotransform[5]
# get target resolution
resolution = target_geotransform[1] # or geotransform[5]
print(str(resolution))
#get target projection
# get target projection
spatial_ref = osr.SpatialReference()
spatial_ref.ImportFromWkt(target_projection)
spatial_ref_projection = spatial_ref.ExportToProj4()
print(spatial_ref_projection)
#gdalwarp call
subprocess.check_output("gdalwarp -dstnodata -32768 -tr "+str(resolution)+" "+str(resolution)+" -r cubicspline -te "+te+" -t_srs '"+spatial_ref_projection+"' "+psrtm+" "+pout, stderr=subprocess.STDOUT, shell=True)
# gdalwarp call
subprocess.check_output(
"gdalwarp -dstnodata -32768 -tr " +
str(resolution) +
" " +
str(resolution) +
" -r cubicspline -te " +
te +
" -t_srs '" +
spatial_ref_projection +
"' " +
psrtm +
" " +
pout,
stderr=subprocess.STDOUT,
shell=True)
def main(argv):
#parse files path
psrtm = argv[1]
pimg = argv[2]
pout = argv[3]
build_dem(psrtm, pimg, pout)
# parse files path
psrtm = argv[1]
pimg = argv[2]
pout = argv[3]
build_dem(psrtm, pimg, pout)
if __name__ == "__main__":
if len(sys.argv) != 4 :
show_help()
else:
main(sys.argv)
if len(sys.argv) != 4:
show_help()
else:
main(sys.argv)
import sys
import subprocess
import subprocess
import glob
import os
import os.path as op
......@@ -12,103 +12,151 @@ import gdalconst
import numpy as np
# run subprocess and write to stdout and stderr
def call_subprocess(process_list):
process = subprocess.Popen(process_list, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
process = subprocess.Popen(
process_list,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE)
out, err = process.communicate()
print out
sys.stderr.write(err)
def get_raster_as_array(raster_file_name):
dataset = gdal.Open(raster_file_name, gdalconst.GA_ReadOnly)
wide = dataset.RasterXSize
high = dataset.RasterYSize
band = dataset.GetRasterBand(1)
array = band.ReadAsArray(0, 0, wide, high)
return array, dataset
dataset = gdal.Open(raster_file_name, gdalconst.GA_ReadOnly)
wide = dataset.RasterXSize
high = dataset.RasterYSize
band = dataset.GetRasterBand(1)
array = band.ReadAsArray(0, 0, wide, high)
return array, dataset
def compute_cloudpercent(image_path):
array_image, dataset_image = get_raster_as_array(image_path)
cloud = np.sum(array_image == 205)
tot_pix = np.sum(array_image != 254)
return (float(cloud)/float(tot_pix))*100
array_image, dataset_image = get_raster_as_array(image_path)
cloud = np.sum(array_image == 205)
tot_pix = np.sum(array_image != 254)
return (float(cloud) / float(tot_pix)) * 100
def compute_snowpercent(image_path):
array_image, dataset_image = get_raster_as_array(image_path)
cloud = np.sum(array_image == 100)
tot_pix = np.sum(array_image != 254)
return (float(cloud)/float(tot_pix))*100
array_image, dataset_image = get_raster_as_array(image_path)
cloud = np.sum(array_image == 100)
tot_pix = np.sum(array_image != 254)
return (float(cloud) / float(tot_pix)) * 100
def format_LIS(snow_detector):
path_img = snow_detector.img
pout = snow_detector.path_tmp
zs = snow_detector.zs
ram = snow_detector.ram
mode = snow_detector.mode
nodata_path = snow_detector.nodata_path
product_id = op.splitext(op.basename(path_img))[0]
path_products = op.join(pout, "LIS_PRODUCTS")
if not op.exists(path_products):
os.makedirs(path_products)
ext = "TIF"
code_snow_all = "_SNOW_ALL"
str_snow_all = product_id+code_snow_all+"."+ext
str_snow_all = str_snow_all.upper()
copyfile(op.join(pout, "snow_all.tif"), op.join(path_products, str_snow_all))
code_seb = "_SEB"
str_seb = product_id+code_seb+"."+ext
str_seb = str_seb.upper()
format_SEB_values(op.join(pout, "final_mask.tif"), nodata_path, ram)
copyfile(op.join(pout, "final_mask.tif"), op.join(path_products, str_seb))
code_seb_vec = "_SEB_VEC"
for f in glob.glob(op.join(pout, "final_mask_vec.*")):
extension = op.splitext(f)[1]
str_seb_vec = product_id+code_seb_vec+extension
str_seb_vec = str_seb_vec.upper()
if extension == ".dbf":
format_SEB_VEC_values(f)
copyfile(f, op.join(path_products, str_seb_vec))
code_compo = "_COMPO"
str_compo = product_id+code_compo+".tif"
str_compo = str_compo.upper()
copyfile(op.join(pout, "composition.tif"), op.join(path_products, str_compo))
#Copy histogram to LIS_PRODUCTS directory
ext = "TXT"
code_histo = "_HISTO"
str_histo = product_id+code_histo+"."+ext
str_histo = str_histo.upper()
copyfile(op.join(pout, "histogram.txt"), op.join(path_products, str_histo))
snow_percent = compute_snowpercent(op.join(pout, "final_mask.tif"))
cloud_percent = compute_cloudpercent(op.join(pout, "final_mask.tif"))
root = etree.Element("Source_Product")
etree.SubElement(root, "PRODUCT_ID").text = product_id
egil = etree.SubElement(root, "Global_Index_List")
etree.SubElement(egil, "QUALITY_INDEX", name='ZS').text = str(zs)
etree.SubElement(egil, "QUALITY_INDEX", name='SnowPercent').text = str(snow_percent)
etree.SubElement(egil, "QUALITY_INDEX", name='CloudPercent').text = str(cloud_percent)
et = etree.ElementTree(root)
et.write(op.join(pout, "metadata.xml"), pretty_print = True)
code_metadata = "_METADATA"
str_metadata = product_id+code_metadata+".xml"
str_metadata = str_metadata.upper()
copyfile(op.join(pout, "metadata.xml"), op.join(path_products, str_metadata))
path_img = snow_detector.img
pout = snow_detector.path_tmp
zs = snow_detector.zs
ram = snow_detector.ram
mode = snow_detector.mode
nodata_path = snow_detector.nodata_path
product_id = op.splitext(op.basename(path_img))[0]
path_products = op.join(pout, "LIS_PRODUCTS")
if not op.exists(path_products):
os.makedirs(path_products)
ext = "TIF"
code_snow_all = "_SNOW_ALL"
str_snow_all = product_id + code_snow_all + "." + ext
str_snow_all = str_snow_all.upper()
copyfile(
op.join(
pout, "snow_all.tif"), op.join(
path_products, str_snow_all))
code_seb = "_SEB"
str_seb = product_id + code_seb + "." + ext
str_seb = str_seb.upper()
format_SEB_values(op.join(pout, "final_mask.tif"), nodata_path, ram)
copyfile(op.join(pout, "final_mask.tif"), op.join(path_products, str_seb))
code_seb_vec = "_SEB_VEC"
for f in glob.glob(op.join(pout, "final_mask_vec.*")):
extension = op.splitext(f)[1]
str_seb_vec = product_id + code_seb_vec + extension
str_seb_vec = str_seb_vec.upper()
if extension == ".dbf":
format_SEB_VEC_values(f)
copyfile(f, op.join(path_products, str_seb_vec))
code_compo = "_COMPO"
str_compo = product_id + code_compo + ".tif"
str_compo = str_compo.upper()
copyfile(
op.join(
pout, "composition.tif"), op.join(
path_products, str_compo))
# Copy histogram to LIS_PRODUCTS directory
ext = "TXT"
code_histo = "_HISTO"
str_histo = product_id + code_histo + "." + ext
str_histo = str_histo.upper()
copyfile(op.join(pout, "histogram.txt"), op.join(path_products, str_histo))
snow_percent = compute_snowpercent(op.join(pout, "final_mask.tif"))
cloud_percent = compute_cloudpercent(op.join(pout, "final_mask.tif"))
root = etree.Element("Source_Product")
etree.SubElement(root, "PRODUCT_ID").text = product_id
egil = etree.SubElement(root, "Global_Index_List")
etree.SubElement(egil, "QUALITY_INDEX", name='ZS').text = str(zs)
etree.SubElement(
egil,
"QUALITY_INDEX",
name='SnowPercent').text = str(snow_percent)
etree.SubElement(
egil,
"QUALITY_INDEX",
name='CloudPercent').text = str(cloud_percent)
et = etree.ElementTree(root)
et.write(op.join(pout, "metadata.xml"), pretty_print=True)
code_metadata = "_METADATA"
str_metadata = product_id + code_metadata + ".xml"
str_metadata = str_metadata.upper()
copyfile(
op.join(
pout, "metadata.xml"), op.join(
path_products, str_metadata))
def format_SEB_values(path, nodata_path, ram):
call_subprocess(["otbcli_BandMath", "-il", path, "-out", path, "uint8", "-ram",str(ram), "-exp", "(im1b1==1)?100:(im1b1==2)?205:0"])
call_subprocess(["otbcli_BandMath", "-il", path, nodata_path, "-out", path, "uint8", "-ram" , str(ram), "-exp", "im2b1==1?254:im1b1"])
call_subprocess(["otbcli_BandMath",
"-il",
path,
"-out",
path,
"uint8",
"-ram",
str(ram),
"-exp",
"(im1b1==1)?100:(im1b1==2)?205:0"])
call_subprocess(["otbcli_BandMath",
"-il",
path,
nodata_path,
"-out",
path,
"uint8",
"-ram",
str(ram),
"-exp",
"im2b1==1?254:im1b1"])
def format_SEB_VEC_values(path):
table = op.splitext(op.basename(path))[0]
call_subprocess(["ogrinfo", path, "-sql", "ALTER TABLE "+table+" ADD COLUMN type varchar(15)"])
call_subprocess(["ogrinfo", path, "-dialect", "SQLite", "-sql", "UPDATE '"+table+"' SET DN=100, type='snow' WHERE DN=1"])
call_subprocess(["ogrinfo", path, "-dialect", "SQLite", "-sql", "UPDATE '"+table+"' SET DN=205, type='cloud' WHERE DN=2"])
call_subprocess(["ogrinfo", path, "-dialect", "SQLite", "-sql", "UPDATE '"+table+"' SET DN=254, type='no data' WHERE DN != 100 AND DN != 205"])
table = op.splitext(op.basename(path))[0]
call_subprocess(["ogrinfo", path, "-sql", "ALTER TABLE " +
table + " ADD COLUMN type varchar(15)"])
call_subprocess(["ogrinfo", path, "-dialect", "SQLite", "-sql",
"UPDATE '" + table + "' SET DN=100, type='snow' WHERE DN=1"])
call_subprocess(["ogrinfo", path, "-dialect", "SQLite", "-sql",
"UPDATE '" + table + "' SET DN=205, type='cloud' WHERE DN=2"])
call_subprocess(["ogrinfo", path, "-dialect", "SQLite", "-sql", "UPDATE '" +
table + "' SET DN=254, type='no data' WHERE DN != 100 AND DN != 205"])
This diff is collapsed.
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