Add smoothing option using shapely

LiCSBAS.yml

@@ -11,4 +11,5 @@ dependencies:
  - numpy
  - psutil
  - requests
+ - shapely
  - statsmodels

LiCSBAS_requirements.txt

@@ -7,4 +7,5 @@ matplotlib
 numpy
 psutil
 requests
+shapely
 statsmodels

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,19 +80,20 @@ 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
 
  #%% 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)