From d0ba1c6ac30ab5c3da74b877e4771387d098b7df Mon Sep 17 00:00:00 2001 From: Yu Morishita Date: Fri, 26 Feb 2021 11:17:33 +0900 Subject: [PATCH] Add smoothing option using shapely --- LiCSBAS.yml | 1 + LiCSBAS_requirements.txt | 1 + bin/LiCSBAS_contour.py | 81 +++++++++++++++++++++++++++++++--------- 3 files changed, 65 insertions(+), 18 deletions(-) diff --git a/LiCSBAS.yml b/LiCSBAS.yml index 0bb6943..796d28c 100644 --- a/LiCSBAS.yml +++ b/LiCSBAS.yml @@ -11,4 +11,5 @@ dependencies: - numpy - psutil - requests + - shapely - statsmodels diff --git a/LiCSBAS_requirements.txt b/LiCSBAS_requirements.txt index c0f8176..e658f22 100644 --- a/LiCSBAS_requirements.txt +++ b/LiCSBAS_requirements.txt @@ -7,4 +7,5 @@ matplotlib numpy psutil requests +shapely statsmodels diff --git a/bin/LiCSBAS_contour.py b/bin/LiCSBAS_contour.py index b8b07e7..33eb49a 100755 --- a/bin/LiCSBAS_contour.py +++ b/bin/LiCSBAS_contour.py @@ -1,33 +1,48 @@ #!/usr/bin/env python3 """ -v1.0 20200408 Yu Morishita, GSI +v1.1 20210226 Yu Morishita, GSI ======== Overview ======== -This script draws contours from a GeoTIFF file and output a GeoJSON file (with GSImaps style). +This script draws contours from a GeoTIFF file and output a GeoJSON file (with GSI Maps style). ===== Usage ===== -LiCSBAS_contour.py -i geotiff -c cont_int [-q cut_nodes] [-o contfile] [-a attrib] [--nodata float] [--no_zero] [--color_n colorcode] [--color_p colorcode] [--color_0 colorcode] [--width float] [--opacity float] +LiCSBAS_contour.py -i geotiff -c cont_int [-q cut_nodes] [-o contfile] + [-a attrib] [-s smoothing_length] [--nodata float] [--no_zero] + [--color_n colorcode] [--color_p colorcode] [--color_0 colorcode] + [--width float] [--opacity float] -i Input GeoTIFF file -c Contour interval -q Do not draw contours with less nodes than this number (Default: 10) -o Output contour GeoJSON file (Default: [geotiff%.tif].cont.geojson) - -a Name for the attribute (good to include unit) (Default: geotiff file name) + -a Name for the attribute (good to include unit) + (Default: geotiff file name) + -s Smoothing length in km (Default: 0; no smoothing) + Note that the shapely module must be installed to use this option. --nodata Nodata value (Default: nan) --no_zero Do not draw contours with 0 - --color_[n|p|0] Color code of contours with negative, positive, 0 values. - (e.g., --color_n "#0000ff" --color_p "#ff0000", blue for negative and red for positive) - (Default: "#000000" (black) for all) + --color_0 Color code of contours with 0 values (Default: #000000 (black)) + --color_n Color code of contours with negative values. + (Default: #0000ff (blue)) + --color_p Color code of contours with positive values. + (Default: #ff0000 (red)) --width Width of contour lines (Default: 2) --opacity Opacity of contour lines (Default: 0.5) +Note: + - color_[n|p|0], width and opacity have an effect in GSI Maps, not in QGIS. + - Rocommend reducing n_node to < 10000 by -q and -s options for GSIMaps. + """ #%% Change log ''' +v1.1 20210226 Yu Morishita, GSI + - Add -s option + - Change default color v1.0 20200408 Yu Morishita, GSI - Original implementationf ''' @@ -49,13 +64,13 @@ def __init__(self, msg): #%% Main def main(argv=None): - + #%% Check argv if argv == None: argv = sys.argv - + start = time.time() - ver=1.0; date=20200408; author="Y. Morishita" + ver=1.1; date=20210226; author="Y. Morishita" print("\n{} ver{} {} {}".format(os.path.basename(argv[0]), ver, date, author), flush=True) print("{} {}".format(os.path.basename(argv[0]), ' '.join(argv[1:])), flush=True) @@ -65,11 +80,12 @@ def main(argv=None): node_thre = 10 outfile = [] attrib = [] + slen_km = 0 nodata = np.nan no_zero_flag = False - color_n = '#000000' - color_p = '#000000' + color_n = '#0000ff' + color_p = '#ff0000' color_0 = '#000000' opacity = 0.5 width = 2 @@ -77,7 +93,7 @@ def main(argv=None): #%% Read options try: try: - opts, args = getopt.getopt(argv[1:], "hi:o:c:q:a:", ["help", "nodata=", "no_zero", "color_n=", "color_p=", "color_0=", "opacity=", "width="]) + opts, args = getopt.getopt(argv[1:], "hi:o:c:q:a:s:", ["help", "nodata=", "no_zero", "color_n=", "color_p=", "color_0=", "opacity=", "width="]) except getopt.error as msg: raise Usage(msg) for o, a in opts: @@ -94,6 +110,8 @@ def main(argv=None): node_thre = int(a) elif o == '-a': attrib = a + elif o == '-s': + slen_km = float(a) elif o == '--nodata': nodata = float(a) elif o == '--no_zero': @@ -128,11 +146,22 @@ def main(argv=None): outfile = infile.replace('.tif', '.cont.geojson') if not attrib: attrib = infile - + if slen_km != 0: + try: + from shapely.geometry import LineString + slen_deg = slen_km*360/40000 #approx + print("\nSmoothing by ~{} km".format(slen_km)) + except ModuleNotFoundError: + print("\nWarning: No module named 'shapely' found") + print(" -s option needs shapely module.") + print(" If you want to use -s, install shapely.") + print(" Deactivate -s option.\n") + slen_km = 0 + #%% Make contour geojson call = ["gdal_contour", "-snodata", str(nodata), "-a", attrib, "-i", str(cint), "-f", "GeoJSON", infile, outfile ] - + p = subp.Popen(call, stdout = subp.PIPE, stderr = subp.STDOUT) for line in iter(p.stdout.readline, b''): print(line.rstrip().decode("utf8")) @@ -143,6 +172,9 @@ def main(argv=None): json_dict = json.load(f) features_list = json_dict['features'] ## list n_feature_in = len(features_list) + n_node_in = np.array([len(feature['geometry']['coordinates']) for + feature in features_list]).sum() + #%% Prepare output features_out_list = [] @@ -151,6 +183,14 @@ def main(argv=None): if no_zero_flag and feature['properties'][attrib] == 0: continue + ## Smoothing + if slen_km != 0: + line = LineString([tuple(this) for this in + feature['geometry']['coordinates']]) + line2 = line.simplify(slen_deg) + feature['geometry']['coordinates'] = \ + [list(i) for i in list(line2.coords)] + ## Remove lines with small n_node n_node = len(feature['geometry']['coordinates']) if n_node <= node_thre: @@ -167,7 +207,7 @@ def main(argv=None): else: feature['properties']['_color'] = color_n - ## Set opacity and weight + ## Set opacity and weight feature['properties']['_opacity'] = opacity feature['properties']['_weight'] = width @@ -176,13 +216,18 @@ def main(argv=None): #%% Output json n_feature_out = len(features_out_list) - print('\nNumber of features: {} -> {}'.format(n_feature_in, n_feature_out)) + n_node_out = np.array([len(feature['geometry']['coordinates']) for + feature in features_out_list]).sum() + print('\nNumber of features: {} -> {}'.format(n_feature_in, + n_feature_out)) + print('Number of nodes: {} -> {}'.format(n_node_in, n_node_out)) if n_feature_out == 0: print('No features remain, not output {}\n'.format(outfile)) os.remove(outfile) else: - jsonout_dict = {'type':json_dict['type'], 'features':features_out_list} + jsonout_dict = {'type':json_dict['type'], + 'features':features_out_list} with open(outfile, 'w') as f: json.dump(jsonout_dict, f, indent=None)