Skip to content

GitLab

Commit 31cab8f7 authored by Rabotin Michael's avatar Rabotin Michael
Browse files

factorized utils method

parent ef4467e9
Hide whitespace changes
Inline Side-by-side
modules/grassUtils.py
View file @ 31cab8f7
...@@ -25,6 +25,7 @@ except ImportError: ...@@ -25,6 +25,7 @@ except ImportError:
#DEVNULL = open('trace_step3', 'wb') #DEVNULL = open('trace_step3', 'wb')
DEVNULL = open(os.devnull, 'wb') DEVNULL = open(os.devnull, 'wb')
import grass.script as grass import grass.script as grass
from grass.script.utils import decode, encode
import subprocess import subprocess
import platform import platform
isWindows = (platform.system() == 'Windows') isWindows = (platform.system() == 'Windows')
...@@ -136,3 +137,46 @@ t-b resol: 1 ...@@ -136,3 +137,46 @@ t-b resol: 1
rc.write('GISDBASE: %s\n' % (os.path.abspath(gisdbase))) rc.write('GISDBASE: %s\n' % (os.path.abspath(gisdbase)))
rc.write('LOCATION_NAME: %s\n' % location) rc.write('LOCATION_NAME: %s\n' % location)
rc.write('MAPSET: %s\n' % mapset) rc.write('MAPSET: %s\n' % mapset)
def reclass(map_in, map_out, size):
p = grass_pipe_command('r.stats', flags='lnNc', input=map_in)
pReclass = grass_feed_command('r.reclass', overwrite=True, input=map_in, output=map_out, rules='-')
for l in p.stdout:
lSpl = decode(l).rstrip(os.linesep).split()
if int(lSpl[1]) <= int(size):
pReclass.stdin.write(encode('%s = %s\n' % (lSpl[0], lSpl[0])))
p.wait()
pReclass.stdin.close()
pReclass.wait()
def count_only(map_in, size):
statLines = decode(grass_read_command('r.stats', quiet=True, flags='lnNc', input=map_in)).rstrip(os.linesep).split(os.linesep)
c = 0
for l in statLines:
lSpl = l.split()
if int(lSpl[1]) <= int(size):
c += 1
return c
#return int(os.popen('r.stats --quiet -lnNc input=%s | awk %s' % (map_in, AWK_cmd_2(size))).readlines()[0])
def getAreasUpTo(map_in, size):
areas = []
statLists = []
statLines = decode(grass_read_command('r.stats', quiet=True, flags='lnNc', input=map_in)).rstrip(os.linesep).split(os.linesep)
for l in statLines:
statLists.append(l.strip().split())
# sort on second column
statLists = sorted(statLists, key=lambda x: int(x[1]))
# awk
intSize = int(size)
old = 0
for sl in statLists:
intSl1 = int(sl[1])
if old != intSl1 and intSl1 <= intSize:
areas.append(intSl1)
old = intSl1
return areas
modules/hrudelin_1_init.py
View file @ 31cab8f7
...@@ -39,133 +39,15 @@ from grass.script.utils import decode, encode ...@@ -39,133 +39,15 @@ from grass.script.utils import decode, encode
from osgeo import ogr from osgeo import ogr
from utils import cutting_raster, get_raster_bounds, get_polygon_bounds
from utils import isfloat, isint, get_coords, is_valid_shp, is_valid_geometry
from utils import is_column_exist, is_column_valid, is_column_value_unique
from utils import are_columns_value_unique
MY_ABS_PATH=os.path.abspath(__file__) MY_ABS_PATH=os.path.abspath(__file__)
MY_DIR=os.path.dirname(MY_ABS_PATH) MY_DIR=os.path.dirname(MY_ABS_PATH)
DATA_PATH=os.path.join(MY_DIR, '..', 'data') DATA_PATH=os.path.join(MY_DIR, '..', 'data')
# cutting the rasters according to user selection
def cutting_raster(raster,raster_out, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
ds = gdal.Open(raster)
ds = gdal.Translate(raster_out, ds, projWin=[xmin_slect, ymax_slect, xmax_slect, ymin_slect])
ds = None
def get_raster_bounds(layer, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
# this is better than sketchy "gdal.Info(ds).split('\n')"
src = gdal.Open(layer)
xmin, xres, xskew, ymax, yskew, yres = src.GetGeoTransform()
xmax = xmin + (src.RasterXSize * xres)
ymin = ymax + (src.RasterYSize * yres)
src = None
return (
max(xmin_slect, xmin),
min(ymax_slect, ymax),
min(xmax_slect, xmax),
max(ymin_slect, ymin)
)
def get_polygon_bounds(shape, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
ds = ogr.Open(shape)
layer = ds.GetLayer()
xmin, xmax, ymin, ymax = layer.GetExtent()
return (
max(xmin_slect, xmin),
min(ymax_slect, ymax),
min(xmax_slect, xmax),
max(ymin_slect, ymin)
)
def isfloat(value):
try:
float(value)
return True
except ValueError:
return False
def isint(value):
try:
int(value)
return True
except ValueError:
return False
def get_coords(parms, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
if not isfloat(parms.get('surface', 'west')) or isfloat(parms.get('surface', 'north')) or isfloat(parms.get('surface', 'east')) or isfloat(parms.get('surface', 'south')):
sys.exit('------------> ERROR : coordinates for surface selection is not valid')
xmin_slect = max(xmin_slect, float(parms.get('surface', 'west')))
ymax_slect = min(ymax_slect, float(parms.get('surface', 'north')))
xmax_slect = min(xmax_slect, float(parms.get('surface', 'east')))
ymin_slect = max(ymin_slect, float(parms.get('surface', 'south')))
print("hello",xmin_slect)
return xmin_slect, ymax_slect, xmax_slect, ymin_slect
def is_valid_shp(shplayer):
if not os.path.isfile(shplayer):
msg="------------> ERROR : Input "+ shplayer+" not found"
sys.exit(msg)
if ogr.Open(shplayer) is None:
msg="------------> ERROR : Input "+ shplayer+" is not a valid shapefile"
sys.exit(msg)
def is_valid_geometry(Name,GeometryType):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
NameLayer_defn=NameLayer.GetLayerDefn()
if ogr.GeometryTypeToName(NameLayer_defn.GetGeomType()) != GeometryType:
msg="------------> ERROR : Input "+ Name+" is not "+ GeometryType
sys.exit(msg)
def is_column_exist(Name,colname):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
NameLayer_defn=NameLayer.GetLayerDefn()
FieldFound=False
for n in range(NameLayer_defn.GetFieldCount()):
fdefn=NameLayer_defn.GetFieldDefn(n)
if fdefn.name == colname:
FieldFound=True
if not FieldFound:
msg="------------> ERROR : column "+ colname+ " not found in Input "+Name
sys.exit(msg)
def is_column_valid(Name,colname):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
listvalue= [ feature.GetField(colname) for feature in NameLayer ]
if None in listvalue or (min(listvalue) == 0):
msg="ERROR : in column "+ colname+ " in Input "+Name +" null or zero value detected"
sys.exit(msg)
def is_column_value_unique(Name,colname):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
listvalue= [ feature.GetField(colname) for feature in NameLayer ]
if(len(set(listvalue)) != len(listvalue)):
msg="ERROR : in column "+ colname+ " in Input "+Name +" non unique value detected"
sys.exit(msg)
def are_columns_value_unique(Name1,colname1,Name2,colname2):
NameDs1=ogr.Open(Name1)
NameLayer1=NameDs1.GetLayer()
listvalue1= [ feature.GetField(colname1) for feature in NameLayer1 ]
NameDs2=ogr.Open(Name2)
NameLayer2=NameDs2.GetLayer()
listvalue2= [ feature.GetField(colname2) for feature in NameLayer2 ]
listvalue1.extend(listvalue2)
if(len(set(listvalue1)) != len(listvalue1)):
msg="ERROR : in column "+ colname1+ " in Input "+Name1 +" and in column "+colname2+" in Input "+Name2+" duplicate value !"
sys.exit(msg)
try: try:
...@@ -309,7 +191,7 @@ def main(parms_file): ...@@ -309,7 +191,7 @@ def main(parms_file):
isinstance(float(irrig_distance_GU),float) isinstance(float(irrig_distance_GU),float)
except: except:
sys.exit('------------> ERROR : irrig_distance_GU PARAMETER is not integer or double') sys.exit('------------> ERROR : irrig_distance_GU PARAMETER is not integer or double')
if irrig_distance_GU<0: if float(irrig_distance_GU)<0:
sys.exit('------------> ERROR : irrig_distance_GU PARAMETER cannot be negative') sys.exit('------------> ERROR : irrig_distance_GU PARAMETER cannot be negative')
......
modules/hrudelin_2_basins.py
View file @ 31cab8f7
...@@ -43,7 +43,9 @@ from osgeo import ogr ...@@ -43,7 +43,9 @@ from osgeo import ogr
import multiprocessing import multiprocessing
from multiprocessing import Pool, cpu_count from multiprocessing import Pool, cpu_count
from utils import isint, write_log
from reach import snapping_points_to_reaches, cut_streams_at_points
from reach import updateAttributeTable, processReachStats
MY_ABS_PATH=os.path.abspath(__file__) MY_ABS_PATH=os.path.abspath(__file__)
MY_DIR=os.path.dirname(MY_ABS_PATH) MY_DIR=os.path.dirname(MY_ABS_PATH)
...@@ -55,395 +57,6 @@ except ImportError: ...@@ -55,395 +57,6 @@ except ImportError:
DEVNULL = open(os.devnull, 'wb') DEVNULL = open(os.devnull, 'wb')
def isint(value):
try:
int(value)
return True
except ValueError:
return False
def write_log(logf,msg, start, end, best_x, best_y,point_code,point_type,point_area):
horiz_shift = best_x - start - end
vertic_shift = best_y - start - end + 1
if horiz_shift != 0 or vertic_shift != 0:
if horiz_shift > 0:
dir_x_shift = 'East'
else:
dir_x_shift = 'West'
if vertic_shift > 0:
dir_y_shift = 'South'
else:
dir_y_shift = 'North'
logf.writerow({
'point code': point_code,
'point type': point_type,
'point area': str(point_area),
'rule': msg,
'horiz dir': dir_x_shift,
'horiz shift': str(abs(horiz_shift)),
'vertic dir': dir_y_shift,
'vertic shift': str(abs(vertic_shift))
})
def snapping_points_to_reaches(parms, directory_out,col_name,col_area,reloc_lyr,reloc_file_name,pointType):
"""
Relocates gauges or dams on the reaches derivated from MNT in step 1 (r.watershed)
"""
basin_min_size = int(parms.get('basin_min_size', 'size'))
surface_ok_1 = int(parms.get('auto_relocation', 'surface_tolerance_1'))
dist_ok_1 = int(parms.get('auto_relocation', 'distance_tolerance_1' ))
if parms.get('auto_relocation', 'surface_tolerance_2') == '':
nb_trials = 1
surface_ok_2 = 0
dist_ok_2 = 0
else:
nb_trials = 2
if not isint(parms.get('auto_relocation', 'surface_tolerance_2')):
sys.exit('------------> ERROR : Surface_tolerance_2 is not integer' )
if not isint(parms.get('auto_relocation', 'distance_tolerance_2')):
sys.exit('------------> ERROR : Distance_tolerance_2 is not integer' )
surface_ok_2 = int(parms.get('auto_relocation', 'surface_tolerance_2'))
dist_ok_2 = int(parms.get('auto_relocation', 'distance_tolerance_2' ))
area_unit = parms.get('auto_relocation', 'area_unit')
# are these the right factors?
# if values is entered in square meters, fact should be 1, no?
if area_unit == '1':
fact = 1000
unit = 'm2'
elif area_unit == '2':
fact = 1000000
unit = 'km2'
else:
print(' =========> ERROR ')
sys.exit(' =========> CHECK area_unit PARAMETER')
logfile_n = os.path.join(directory_out, reloc_file_name)
logfile = open(logfile_n, 'w')
fieldnames = ['point code', 'point type','point area', 'rule', 'horiz shift', 'horiz dir', 'vertic shift', 'vertic dir']
logf = csv.DictWriter(logfile, fieldnames=fieldnames, delimiter=';')
logf.writeheader()
# Read the accumulation raster and get the geo informations
accum_file = os.path.join(directory_out, 'step1_accum.tif')
accum_lyr = gdal.Open(accum_file, GA_ReadOnly)
accum_bd = accum_lyr.GetRasterBand(1)
georef = accum_lyr.GetGeoTransform()
hres = georef[1]
vres = georef[5]
pixel_sz = (hres * abs(vres))
# It gets the number of rows and columns of the accumulation raster file
acc_rast_cols = accum_lyr.RasterXSize
acc_rast_rows = accum_lyr.RasterYSize
# Matrix of euclidian distances from the center (gauge position)
max_dist = max(dist_ok_1, dist_ok_2)
euclid_dist = np.full((max_dist * 2 + 1, max_dist * 2 + 1), 0, dtype=int)
for y in range(0, (max_dist * 2 + 1)):
for x in range(0, (max_dist * 2 + 1)):
euclid_dist[y, x] = int(math.sqrt((y - max_dist)**2 + (x - max_dist)**2))
# Points loop
for point in reloc_lyr:
if not point.IsFieldSet(col_name):
sys.exit('------------> ERROR : col_name not found' )
point_code = point.GetField(col_name)
# Get the drained surface of the point
point_area = point.GetField(col_area)
# ------------------------------------------------------------
# if point_area < 0:
# nb_pixel_area = basin_min_size + 1
# ignore_surface_tolerance = True
# else:
# nb_pixel_area = ((point_area * fact) / pixel_sz) + 0.5
# ignore_surface_tolerance = False
# ------------------------------------------------------------
# here, the corresponding number of cell of the specified drained area is calculated
nb_pixel_area = ((point_area * fact) / pixel_sz) + 0.5
# ------------------------------------------------------------
if nb_pixel_area > basin_min_size:
# Get the coordinates of the point and calculate center of the offset
# I think it relocates the station in the center of the closest pixel
# in the raster
geom = point.GetGeometryRef()
point_x, point_y = geom.GetX(), geom.GetY()
Xreste = point_x % hres
Yreste = point_y % abs(vres)
new_point_x = (point_x - Xreste) + 1
new_point_y = (point_y - Yreste) - 1
# Read the currents of accum raster around the point
# (new_point_x - georef[0]) / hres)
# ==> this gives the column in raster file where the station has been relocated
# int(((new_point_x - georef[0]) / hres) - (dist_max))
# ==> get the difference with the maximum number of column of accepted differences
# if > 0 then, dist_max is < column number of point
# if < 0 then, dist_max is > column number of point
# max(..., 0) ==> there, we just see if point is or not within a dist_max distance from left border???
dist_max = max(dist_ok_1, dist_ok_2)
pix_x = max(int(((new_point_x - georef[0]) / hres) - (dist_max)), 0)
pix_y = max(int(((georef[3] - new_point_y) / abs(vres)) - (dist_max)), 0)
# and don't access beyond the accum raster bounds
nb_cols = min((dist_max * 2) + 1, acc_rast_cols - pix_x)
nb_rows = (min((dist_max * 2) + 1, acc_rast_rows - pix_y))
zone = accum_bd.ReadAsArray(pix_x, pix_y, nb_cols, nb_rows).astype(int)
# in the end, this is just to ensure that we won't look at relocating the gaug outside the raster boudaries
# isn't that right?
# zone contains the accumulation raster data within max_dist area
# Matrix of differences
zone_diff = abs(nb_pixel_area - abs(zone))
# here, the differences between the actual drained area and the specified drained area is calculated
# Prepare for first search
# here accum_toler is the number of difference in pixel that can be tolerated
accum_toler = int(nb_pixel_area * surface_ok_1 / 100)
if dist_ok_1 > dist_ok_2:
start = 0
col_max = min(nb_cols, acc_rast_cols - pix_x)
row_max = min(nb_rows, acc_rast_rows - pix_y)
end = dist_ok_1 - start
else:
start = dist_ok_2 - dist_ok_1
col_max = min(nb_cols - start, acc_rast_cols - pix_x)
row_max = min(nb_rows - start, acc_rast_rows - pix_y)
end = dist_ok_2 - start
msg = 'shifted with first rule'
trial = 1
while trial <= nb_trials:
# best_area is useless... not used anywhere...
best_area = zone_diff.max()
best_euclid_dist = euclid_dist[start, start]
best_x = 0
best_y = 0
better_place_found = False
for y in range(start, row_max):
for x in range(start, col_max):
# ------------------------------------------------------------
# if ignore_surface_tolerance:
# if euclid_dist[y, x] < best_euclid_dist:
# best_euclid_dist = euclid_dist[y, x]
# best_y = y
# best_x = x
# better_place_found = True
# else:
# if zone_diff[y, x] <= accum_toler:
# if euclid_dist[y, x] < best_euclid_dist:
# best_euclid_dist = euclid_dist[y, x]
# best_y = y
# best_x = x
# better_place_found = True
# ------------------------------------------------------------
if zone_diff[y, x] <= accum_toler:
if euclid_dist[y, x] < best_euclid_dist:
# best_area is useless... not used anywhere...
best_area = zone_diff[y, x]
best_euclid_dist = euclid_dist[y, x]
best_y = y
best_x = x
better_place_found = True
# ------------------------------------------------------------
if better_place_found:
new_x = ((pix_x + best_x) * hres) + georef[0] + (hres / 2)
new_y = georef[3] - ((pix_y + best_y) * abs(vres)) - (abs(vres) / 2)
wkt = 'POINT(%f %f)' % (new_x , new_y)
pointWkt = ogr.CreateGeometryFromWkt(wkt)
point.SetGeometryDirectly(pointWkt)
reloc_lyr.SetFeature(point)
trial = 3
write_log(logf,msg, start, end, best_x, best_y,point_code,pointType,point_area)
else:
if trial == 2:
logf.writerow({'point code': point_code, 'point type':pointType,'point area': str(point_area), 'rule': ' no better place found'})
trial = 3
else:
if nb_trials == 1:
logf.writerow({'point code': point_code, 'point type':pointType,'point area': str(point_area), 'rule': ' no better place found'})
trial = 3
else:
# Prepare for second search
msg = 'shifted with second rule'
accum_toler = int(nb_pixel_area * surface_ok_2 / 100)
trial = 2
if dist_ok_2 > dist_ok_1:
start = 0
end = dist_ok_2 - start
else:
start = dist_ok_1 - dist_ok_2
col_max = min(nb_cols - start, acc_rast_cols - pix_x)
row_max = min(nb_rows - start, acc_rast_rows - pix_y )
end = dist_ok_1 - start
else:
pointID = point.GetFID()
reloc_lyr.DeleteFeature(pointID)
logf.writerow({'point code': point_code, 'point type':pointType,'point area': str(point_area), 'rule': ' basin too small ---> DELETED'})
logfile.close
def cut_streams_at_points(directory_out, list_point,stream_in,stream_out):
listReachID=[]
"""
Streams generated at step 1 (by r.watershed) are cutted at the gauges or at the dams
"""
streams_file = os.path.join(directory_out, stream_in)
streams_lyr = gdal.Open(streams_file, GA_ReadOnly)
streams_band = streams_lyr.GetRasterBand(1)
stream_orig = streams_band.ReadAsArray()
mask_file = os.path.join(directory_out, 'step2_mask.tif')
mask_lyr = gdal.Open(mask_file, GA_ReadOnly)
mask_band = mask_lyr.GetRasterBand(1)
mask = mask_band.ReadAsArray()
nodata = mask_band.GetNoDataValue()
mask[mask == nodata] = 0
stream_orig = stream_orig * mask
streams_cols = streams_lyr.RasterXSize
streams_rows = streams_lyr.RasterYSize
georef = streams_lyr.GetGeoTransform()
hres = georef[1]
vres = georef[5]
streams_new_file = os.path.join(directory_out, stream_out)
driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create(streams_new_file, streams_cols, streams_rows, 1, GDT_Int32)
dst_ds.SetGeoTransform((georef[0], hres, 0, georef[3], 0, vres))
dst_ds.SetProjection(streams_lyr.GetProjection())
drain_file = os.path.join(directory_out, 'step1_drain.tif')
drain_lyr = gdal.Open(drain_file, GA_ReadOnly)
drain = drain_lyr.GetRasterBand(1).ReadAsArray()
drain_dir = [[-1,1],[-1,0],[-1,-1],[0,-1],[1,-1],[1,0],[1,1],[0,1]]
# Give possibility of more than 2 segments in a stream
stream_wk = np.zeros((streams_rows, streams_cols), dtype=('int64'))
nodata = streams_band.GetNoDataValue()
for y in range (0, streams_rows - 1):
for x in range (0, streams_cols - 1):
if stream_orig[y, x] != nodata:
stream_wk[y, x] = stream_orig[y, x] * 100
for point in list_point:
current_x = int(((point[0] - georef[0]) / hres))
current_y = int(((georef[3] - point[1]) / abs(vres)))
current_stream = stream_orig[current_y, current_x]
if current_stream == nodata:
print('---------- gauge outside of stream network: ', point)
end_of_stream = True
else:
stream_wk[current_y, current_x] += 1
end_of_stream = False
i = 0
d = abs(drain[current_y, current_x])
while not end_of_stream:
if d == 8:
d = 0
y = current_y + drain_dir[d][0]
x = current_x + drain_dir[d][1]
d2 = abs(drain[y, x]) - 1
drain_into_y = y + drain_dir[d2][0]
drain_into_x = x + drain_dir[d2][1]
if (stream_orig[y, x] == current_stream) and (drain_into_y == current_y) and (drain_into_x == current_x):
stream_wk[y, x] += 1
current_y = y
current_x = x
i = 0
d = abs(drain[y, x])
else:
i += 1
if i == 7:
end_of_stream = True
else:
d += 1
dst_ds.GetRasterBand(1).WriteArray(stream_wk)
dst_ds = None
for point in list_point:
current_x = int(((point[0] - georef[0]) / hres))
current_y = int(((georef[3] - point[1]) / abs(vres)))
current_stream = stream_wk[current_y, current_x]
tuple_ReachID=(point[2],current_stream)
listReachID.append(tuple_ReachID)
return listReachID
def updateAttributeTable(directory_out,list_point,shapefile,colname):
point_file = os.path.join(directory_out, shapefile)
point_in = ogr.Open(point_file, 1)
point_lyr = point_in.GetLayer()
point_lyr.CreateField(ogr.FieldDefn('ReachID', ogr.OFTInteger))
point_features = point_lyr.GetNextFeature()
while point_features:
pointID=point_features.GetField(colname)
reachID=0
for point in list_point: