[GRASS-SVN] r69067 - in grass-addons/grass7/imagery: . i.image.bathymetry
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Aug 3 05:28:19 PDT 2016
Author: vinayocu
Date: 2016-08-03 05:28:19 -0700 (Wed, 03 Aug 2016)
New Revision: 69067
Added:
grass-addons/grass7/imagery/i.image.bathymetry/
grass-addons/grass7/imagery/i.image.bathymetry/Makefile
grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html
grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py
Log:
SDB
Added: grass-addons/grass7/imagery/i.image.bathymetry/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.image.bathymetry/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.image.bathymetry/Makefile 2016-08-03 12:28:19 UTC (rev 69067)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.image.bathymetry
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/imagery/i.image.bathymetry/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html
===================================================================
--- grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html (rev 0)
+++ grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html 2016-08-03 12:28:19 UTC (rev 69067)
@@ -0,0 +1,32 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.image.bathymetry</em> is used to estimate near-shore depth from optical satellite images. Module estimates bathymetry over near-shore region using limited reference depth points. The maximum depth can be estimated by the module is depending up on many factors such as quality of the water and suspended materials etc. (Lyzenga et al., 2006, Kanno and Tanaka, 2012). Our experiments with several multi-spectral optical images indicate that the depth estimates are reliable for when water column is below 20 meter.
+<br><br>
+Delineation of land and water area are based on combining the result of NDVI and band ratio. NDVI has used to delineate water from land, band ratio between green band and infrared band used to separate the delineated water from clouds, ice etc. Atmospheric and water column correction applied according to the Lyzenga et al., 2006. Corrected spectral bands will be used for weighted multiple regression to estimate depth. R library <em>GWmodel</em> has been used to compute the Geographically Weighted Regression used for depth estimation.
+<h2>NOTES</h2>
+The input image must include deep water pixels (far away from the coast) which are used to assist water surface and water column correction. Sparse depth points extracted from hydrographic charts or depth pints derived from LiDAR survey or derived from Sonar survey can be used as reference depth for calibration.
+<br><br>
+ The tide height at the time of reference depth collection and satellite imagery capture should be normalized if it is not. An option is available in the module to provide tide hieght at the tide of image captured and the module will correct the reference depth accordingly. This option asuumes that the reference depth given is corrected zero tide height. The tide lower than zero can be added as negative value.
+<br><br>
+The <em>GWmodel</em> adaptive GWR model is memory intensive and cannot be used to process large images. For large images, the estimation is carried out by using non-adaptive GWR implemented in <em>r.gwr</em> module in GRASS GIS. R > 3.1 should be installed to run <em>GWmodel</em> in order to proccess adaptive GWR model for better depth estimation. Default gaussian kernel will be used to estimate geographically weighted regression coefficients.The flag 'b' can be used to change the kernel function gaussian to bi-square.
+<br><br>
+<h2>EXAMPLES</h2>
+In <em>i.image.bathymetry</em> green band, red band, near-infrared band, band for correction and calibration depth points are mandatory input. Additional bands available in the visible wavelength can be used for better depth estimation as optional input. Short Wave Infrared (SWIR) band is suggested to use as "band_for_correction" if it is available (for e.g. satellite images like Landsat-7, Landsat-8 and Sentinel-2).An example of depth estimation using Sentinel-2 (MSI) image is shown below
+<br><br>
+i.image.bathymetry band1_blue='B2' band2_green='B3' band3_red='B4' band5_nir='B8' band_for_correction='B11' calibration_points='Calibration_points' depth_estimate='output'
+<br><br>
+If SWIR band is not available near-infrared band can be used as "band_for_correction" (for e.g. satellite images like RapidEye and ALOS AVINIR-2). An example of depth estimation using RapidEye image is shown below image is shown below
+<br><br>
+i.image.bathymetry band1_blue='B1' band2_green='B2' band3_red='B3' Additional_band1='B4' band5_nir='B5' band_for_correction='B5' calibration_points='Calibration_points' depth_estimate='output'
+<h2>AUTHORS</h2>
+Vinayaraj Poliyapram (email:<em>vinay223333 at gmail.com</em>), Luca Delucchi and Venkatesh Raghavan
+<h2>SEE ALSO</h2>
+<em>r.gwr</em> and <em>r.regression.multi</em>
+<h2>REFERENCES</h2>
+<ul>
+<li>Lyzenga, D. R., Malinas, N. R. and Tanis, F. J., 2006, Multispectral Bathymetry using a Simple Physically Based Algorithm. IEEE Trans. Geosci. Remote Sensing, 44(8): 2251-2259.
+
+<li>Su, H., Liu, H., Lei, W., Philipi, M., Heyman, W., and Beck, A., 2013, Geographically Adaptive Inversion Model for Improving Bathymetric Retrieval from Multispectral satellite Imagery. IEEE Transaction on Geosciences and Remote Sensing, 52(1) : 465-476, Accessed January 2013, doi:10.1109/TGRS.2013.2241772.
+
+<li>Vinayaraj, P., Raghavan, V., Masumoto, S. and Glejin, J., 2015, Comparative Evaluation and Refinement of Algorithm for Water Depth Estimation using Medium Resolution Remote Sensing Data. International Journal of geoinformatics, 11(3):17-29
+
Property changes on: grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py
===================================================================
--- grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py (rev 0)
+++ grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py 2016-08-03 12:28:19 UTC (rev 69067)
@@ -0,0 +1,316 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: i.image.bathymetry
+# AUTHOR(S): Vinayaraj Poliyapram <vinay223333 at gmail.com> and Luca Delulucchi
+#
+# PURPOSE: Script for estimating bathymetry from optical satellite images
+# COPYRIGHT: (C) Vinayaraj Poliyapram and by the GRASS Development Team
+#
+# This program is free software under the GNU General
+# Public License (>=v2). Read the file COPYING that
+# comes with GRASS for details.
+#
+#############################################################################
+
+#%module
+#% description: Satellite Derived Bathymetry (SDB) from multispectral images
+#% keyword: Bathymetry
+#% keyword: Satellite
+#%end
+#%option G_OPT_R_INPUT
+#% key: blue_band
+#%required: no
+#%end
+#%option G_OPT_R_INPUT
+#% key: green_band
+#%required: yes
+#%end
+#%option G_OPT_R_INPUT
+#% key: red_band
+#%required: yes
+#%end
+#%option G_OPT_R_INPUT
+#% key: nir_band
+#%required: yes
+#%end
+#%option G_OPT_R_INPUT
+#% key: band_for_correction
+#%required: yes
+#%end
+#%option G_OPT_V_INPUT
+#% key: calibration_points
+#%required: yes
+#%end
+#%option G_OPT_R_INPUT
+#% key: additional_band1
+#%required: no
+#%end
+#%option G_OPT_R_INPUT
+#% key: additional_band2
+#%required: no
+#%end
+#%option G_OPT_R_INPUT
+#% key: additional_band3
+#%required: no
+#%end
+#%option G_OPT_R_INPUT
+#% key: additional_band4
+#%required: no
+#%end
+#%option G_OPT_R_OUTPUT
+#% key: depth_estimate
+#%required: yes
+#%end
+#%option
+#% key: tide_height
+#%type: double
+#%multiple: no
+#%required: no
+#%description: Tide correction to the time of satellite image capture
+#%end
+#%flag
+#%key: b
+#% description: select kernel function as bi-square
+#%end
+
+import atexit
+import grass.script as g
+from grass.pygrass.raster import RasterRow
+from grass.pygrass.modules import Module
+import subprocess
+import os
+crctd_lst = []
+
+
+def main():
+ options, flags = g.parser()
+ Blue = options['blue_band']
+ Green = options['green_band']
+ Red = options['red_band']
+ NIR = options['nir_band']
+ SWIR = options['band_for_correction']
+ Calibration_points = options['calibration_points']
+ Additional_band1 = options['additional_band1']
+ Additional_band2 = options['additional_band2']
+ Additional_band3 = options['additional_band3']
+ Additional_band4 = options['additional_band4']
+ bathymetry = options['depth_estimate']
+ tide_height = options['tide_height']
+ bisquare = flags['b']
+ g.run_command('g.region', raster=Green)
+ g.run_command('v.to.rast', input=Calibration_points, type='point',
+ use='attr', attribute_column='value',
+ output='tmp_Calibration_points')
+ if tide_height:
+ cal = g.parse_command('r.univar', map='tmp_Calibration_points',
+ flags='g')
+ if float(cal['min']) >= 0:
+ t = float(tide_height)
+ g.mapcalc(exp="{d}=({d}+float({t}))".format(d='tmp_Calibration_points',
+ t=t),
+ overwrite=True)
+ if float(cal['min']) < 0:
+ t = float(tide_height) * -1
+ g.mapcalc(exp="{d}=({d}+float({t}))".format(d='tmp_Calibration_points',
+ t=t),
+ overwrite=True)
+ g.mapcalc(exp="{tmp_ratio}=({Green}/{SWIR})".format(tmp_ratio='tmp_ratio',
+ Green=Green, SWIR=SWIR))
+ g.mapcalc(exp="{tmp_NDVI}=float({NIR}-{Red})/float({NIR}+{Red})"
+ .format(tmp_NDVI='tmp_NDVI', NIR=NIR, Red=Red))
+ g.mapcalc(exp="{tmp_water}=if({tmp_ratio} < 1, null(), if({tmp_NDVI} <"
+ "0, {tmp_ratio}, null()))".format(tmp_NDVI='tmp_NDVI',
+ tmp_water='tmp_water',
+ tmp_ratio='tmp_ratio'))
+ g.run_command('r.mask', raster='tmp_water', overwrite=True)
+ li = [Green, Additional_band1, Additional_band2, Additional_band3,
+ Additional_band4, Blue, Red]
+
+ for i in li:
+ j, sep, tail = i.partition('@')
+ tmp_ = RasterRow(str(i))
+ if tmp_.exist() is False:
+ continue
+ g.run_command('g.region', vector=Calibration_points)
+ # Avoid zero values
+ g.mapcalc(exp="{tmp_b}=if({x}>1, {x},null())".format(tmp_b='tmp_b',
+ x=str(i)), overwrite=True)
+
+ tmp_AOI = g.parse_command('r.univar', map='tmp_b', flags='g')
+ tmp_AOI_min = float(tmp_AOI['min'])
+ g.run_command('g.region', raster=Green)
+ try:
+ g.mapcalc(exp="{tmp_deep}=if({tmp_band}<{band_min}, {tmp_band},"
+ "null())".format(tmp_deep='tmp_deep',
+ band_min=tmp_AOI_min,
+ tmp_band=str(i)),
+ overwrite=True)
+ g.run_command('r.mask', raster='tmp_deep', overwrite=True)
+ tmp_coe = g.parse_command('r.regression.line', mapx=SWIR,
+ mapy=str(i), flags='g')
+ g.run_command('r.mask', raster='tmp_water', overwrite=True)
+ g.run_command('g.region', vector=Calibration_points)
+ g.mapcalc(exp="{tmp_crctd}=log({tmp_band}-{a}-{b}*{SWIR})"
+ .format(tmp_crctd='tmp_crctd' + str(j),
+ tmp_band=str(i), a=float(tmp_coe['a']),
+ b=float(tmp_coe['b']), SWIR=SWIR),
+ overwrite=True)
+ except:
+ g.message("Cannot find deep water pixels")
+ g.run_command('g.remove', type='raster', name='MASK', flags='f')
+ g.run_command('g.region', vector=Calibration_points)
+ g.mapcalc("{tmp_crctd} = log({tmp_band}-{SWIR})"
+ .format(tmp_crctd='tmp_crctd' + str(j), tmp_band=str(i),
+ SWIR=SWIR),
+ overwrite=True)
+ crctd_lst.append('tmp_crctd' + str(j))
+ # For R-GWmodel
+ r_file = open('crctd.R', 'w')
+ libs = ['GWmodel', 'data.table', 'rgrass7', 'rgdal', 'raster']
+ for i in libs:
+ install = 'if(!is.element("%s", installed.packages()[,1])){\n' % i
+ install += "cat('\\n\\nInstalling %s package from CRAN\n')\n" % i
+ install += "if(!file.exists(Sys.getenv('R_LIBS_USER'))){\n"
+ install += "dir.create(Sys.getenv('R_LIBS_USER'), recursive=TRUE)\n"
+ install += ".libPaths(Sys.getenv('R_LIBS_USER'))}\n"
+ install += 'install.packages("%s", repos="http://cran.us.r-' \
+ 'project.org")}\n' % i
+ r_file.write(install)
+ libraries = 'library(%s)\n' % i
+ r_file.write(libraries)
+ Green_new, sep, tail = Green.partition('@')
+ r_file.write('grass_file = readRAST("tmp_crctd%s")\n' % Green_new)
+ r_file.write('raster_file = raster(grass_file)\n')
+ frame_file = 'pred = as.data.frame(raster_file,na.rm = TRUE,xy = TRUE)\n'
+ r_file.write(frame_file)
+ for i in li:
+ j, sep, tail = i.partition('@')
+ Green_new, sep, tail = Green.partition('@')
+ tmp_ = RasterRow(str(i))
+ if tmp_.exist() is False:
+ continue
+ r_file.write('grass_file = readRAST("tmp_crctd%s")\n' % j)
+ r_file.write('raster_file = raster(grass_file)\n')
+ r_file.write('frame_pred%s = as.data.frame(raster_file, na.rm = TRUE,'
+ 'xy = TRUE)\n' % j)
+ pred_file = 'frame_pred_green=data.frame(frame_pred%s)\n' % Green_new
+ pred_file += 'pred=merge(pred, frame_pred%s)\n' % j
+ r_file.write(pred_file)
+ # For reference_file repeat with MASK
+ g.run_command('r.mask', raster='tmp_Calibration_points',
+ overwrite=True)
+ r_file.write('grass_file=readRAST("%s")\n' % 'tmp_Calibration_points')
+ r_file.write('raster_file = raster(grass_file)\n')
+ frame_file = 'calib = as.data.frame(raster_file,na.rm = TRUE ,' \
+ 'xy = TRUE)\n'
+ r_file.write(frame_file)
+ for i in li:
+ j, sep, tail = i.partition('@')
+ tmp_ = RasterRow(str(i))
+ if tmp_.exist() is False:
+ continue
+ r_file.write('grass_file = readRAST("tmp_crctd%s")\n' % j)
+ r_file.write('raster_file = raster(grass_file)\n')
+ r_file.write('frame_ref%s = as.data.frame(raster_file,na.rm = TRUE,' \
+ 'xy = TRUE)\n' % j)
+ ref_file = 'calib = merge(calib, frame_ref%s)\n' % j
+ r_file.write(ref_file)
+ g.run_command('g.remove', type='raster', pattern='MASK', flags='f')
+ try:
+ ref_file = 'Rapid_ref.sdf=SpatialPointsDataFrame(calib[,1:2],calib)\n'
+ ref_file += 'Rapid_pred.sdf=SpatialPointsDataFrame(pred[,1:2],' \
+ 'pred)\n'
+ ref_file += 'DM_Rapid_ref.sdf=gw.dist(dp.locat=coordinates' \
+ '(Rapid_ref.sdf))\n'
+ r_file.write(ref_file)
+ l = []
+ # Join the corrected bands in to a string
+ le = len(crctd_lst)
+ for i in crctd_lst:
+ l.append(i)
+ k = '+'.join(l)
+ if bisquare:
+ ref_flag = "cat('\nRunning adaptive GWR using" \
+ "bisquare kernel..\n')\n"
+ ref_flag += 'BW_Rapid_ref.sdf=bw.gwr(tmp_Calibration_points~%s,' \
+ 'data=Rapid_ref.sdf, kernel="bisquare",' \
+ 'adaptive=TRUE, dMat=DM_Rapid_ref.sdf)\n' % k
+ ref_flag += 'DM_Rapid_pred.sdf=gw.dist(dp.locat=coordinates' \
+ '(Rapid_ref.sdf), rp.locat=coordinates' \
+ '(Rapid_pred.sdf))\n'
+ ref_flag += 'GWR_Rapid_pred.sdf=gwr.predict(tmp_Calibration_poi' \
+ 'nts~%s,data=Rapid_ref.sdf, bw = BW_Rapid_ref.sdf,' \
+ 'predictdata = Rapid_pred.sdf, kernel = "bisquare",' \
+ 'adaptive = TRUE, dMat1 = DM_Rapid_pred.sdf,' \
+ 'dMat2 = DM_Rapid_ref.sdf)\n' % k
+ r_file.write(ref_flag)
+ if not bisquare:
+ ref_fla = "cat('\nRunning adaptive-GWR using gaussian kernel\n')\n"
+ ref_fla += 'BW_Rapid_ref.sdf=bw.gwr(tmp_Calibration_points~%s,' \
+ 'data=Rapid_ref.sdf, kernel="gaussian",' \
+ 'adaptive=TRUE, dMat= DM_Rapid_ref.sdf)\n' % k
+ ref_fla += 'DM_Rapid_pred.sdf=gw.dist(dp.locat=coordinates' \
+ '(Rapid_ref.sdf), rp.locat=coordinates' \
+ '(Rapid_pred.sdf))\n'
+ ref_fla += 'GWR_Rapid_pred.sdf = gwr.predict(tmp_Calibration_poi' \
+ 'nts~%s,data=Rapid_ref.sdf, bw=BW_Rapid_ref.sdf,' \
+ 'predictdata = Rapid_pred.sdf, kernel = "gaussian",' \
+ 'adaptive = TRUE, dMat1 = DM_Rapid_pred.sdf,' \
+ 'dMat2 = DM_Rapid_ref.sdf)\n' % k
+ r_file.write(ref_fla)
+ ref_fil = 'Sp_frame = as.data.frame(GWR_Rapid_pred.sdf$SDF)\n'
+ r_file.write(ref_fil)
+ r_file.write('write.table(Sp_frame, quote=FALSE, sep=",",' \
+ '"prediction.txt")\n')
+ r_file.close()
+ subprocess.check_call(['Rscript', 'crctd.R'], shell=False)
+ g.run_command('r.in.xyz', input='prediction.txt',
+ output='tmp_bathymetry', skip=1, separator=",",
+ x=(int(le) + 5), y=(int(le) + 6), z=(int(le) + 3),
+ overwrite=True)
+ except subprocess.CalledProcessError:
+ g.message("Integer outflow... ")
+ if bisquare:
+ g.message("Running fixed GWR using bisqare kernel...")
+ bw = g.parse_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points',
+ kernel='bisquare', flags='ge')
+ g.run_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points',
+ estimates='tmp_bathymetry', kernel='bisquare',
+ bandwidth=int(bw['estimate']))
+ else:
+ g.message("Running fixed GWR using gaussian kernel...")
+ try:
+ Module('r.gwr')
+ except:
+ g.run_command('g.extension', extension='r.gwr')
+ bw = g.parse_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points', flags='ge')
+ g.run_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points',
+ estimates='tmp_bathymetry',
+ bandwidth=int(bw['estimate']))
+ tmp_rslt_ext = g.parse_command('r.univar', map='tmp_Calibration_points',
+ flags='g')
+ g.mapcalc(exp="{bathymetry}=if({tmp_SDB}>{max_}, null(),"
+ "if({tmp_SDB}<{min_}, null(), {tmp_SDB}))".format
+ (tmp_SDB='tmp_bathymetry', bathymetry=bathymetry,
+ max_=float(tmp_rslt_ext['max']),
+ min_=float(tmp_rslt_ext['min'])))
+
+
+def cleanup():
+ g.run_command('g.remove', type='raster', pattern='MASK', flags='f')
+ g.run_command('g.remove', type='raster', pattern='*tmp*', flags='f')
+ try:
+ os.remove('crctd.R')
+ os.remove('prediction.txt')
+ except OSError:
+ pass
+
+if __name__ == '__main__':
+ atexit.register(cleanup)
+ main()
More information about the grass-commit
mailing list