[GRASS-SVN] r70144 - grass-addons/grass7/imagery/i.image.bathymetry
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 27 05:17:37 PST 2016
Author: vinayocu
Date: 2016-12-27 05:17:37 -0800 (Tue, 27 Dec 2016)
New Revision: 70144
Modified:
grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html
grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py
Log:
added new flag
Modified: grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html
===================================================================
--- grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html 2016-12-27 12:38:16 UTC (rev 70143)
+++ grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.html 2016-12-27 13:17:37 UTC (rev 70144)
@@ -1,111 +1,31 @@
<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.
-
-<p>
-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.
-
-
+<em>i.image.bathymetry</em> is used to estimate Satellite-Derived Bathymetry (SDB). 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.
-
-<p>
-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
-
-<p>
-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>-b</b> can be
-used to change the kernel function gaussian to bi-square.
-
-
+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, where depth value is stored in column named 'Z'.
-
-<p>
-<div class="code"><pre>
-i.image.bathymetry blue_band='B2' green_band='B3' red_band='B4' nir_band='B8' \
- band_for_correction='B11' calibration_points='Calibration_points'
- calibration_column='Z' depth_estimate='output'
-</pre></div>
-
-<p>
-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, where depth value is stored in column named
-'value'.
-
-<p>
-<div class="code"><pre>
-i.image.bathymetry blue_band='B1' green_band='B2' red_band='B3' \
- additional_band1='B4' nir_band='B5' band_for_correction='B5' \
- calibration_points='Calibration_points' calibration_column='value' \
- depth_estimate='output'
-</pre></div>
-
-
+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, where depth value is stored in column named 'Z'
+<br><br>
+i.image.bathymetry blue_band='B2' green_band='B3' red_band='B4' nir_band='B8' band_for_correction='B11' calibration_points='Calibration_points' calibration_column='Z' 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, where depth value is stored in column named 'value'.
+<br><br>
+i.image.bathymetry blue_band='B1' green_band='B2' red_band='B3' Additional_band1='B4' nir_band='B5' band_for_correction='B5' calibration_points='Calibration_points' calibration_column='value' depth_estimate='output'
<h2>AUTHORS</h2>
-Vinayaraj Poliyapram (email:<em>vinay223333 at gmail.com</em>), Luca Delucchi
-and Venkatesh Raghavan
-
+Vinayaraj Poliyapram (email:<em>vinay223333 at gmail.com</em>), Luca Delucchi and Venkatesh Raghavan
<h2>SEE ALSO</h2>
-
-<em>
-<a href="https://grass.osgeo.org/grass70/manuals/addons/r.gwr.html">r.gwr</a>,
-<a href="r.regression.multi.html">r.regression.multi</a>
-</em>
-
+<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>Vinayaraj, P., Raghavan, V. and Masumoto, S. (2016) Satellite derived bathymetry using adaptive-geographically weighted regression model, Marine Geodesy, 39(6), pp.458-478
-<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>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
-</ul>
+
Modified: grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py
===================================================================
--- grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py 2016-12-27 12:38:16 UTC (rev 70143)
+++ grass-addons/grass7/imagery/i.image.bathymetry/i.image.bathymetry.py 2016-12-27 13:17:37 UTC (rev 70144)
@@ -15,71 +15,74 @@
#############################################################################
#%module
-#% description: Satellite Derived Bathymetry (SDB) from multispectral images.
-#% keywords: imagery
-#% keywords: bathymetry
-#% keywords: satellite
+#% description: Satellite Derived Bathymetry (SDB) from multispectral images
+#% keyword: Bathymetry
+#% keyword: Satellite
#%end
#%option G_OPT_R_INPUT
#% key: blue_band
-#% required: no
+#%required: no
#%end
#%option G_OPT_R_INPUT
#% key: green_band
-#% required: yes
+#%required: yes
#%end
#%option G_OPT_R_INPUT
#% key: red_band
-#% required: yes
+#%required: yes
#%end
#%option G_OPT_R_INPUT
#% key: nir_band
-#% required: yes
+#%required: yes
#%end
#%option G_OPT_R_INPUT
#% key: band_for_correction
-#% required: yes
+#%required: yes
#%end
#%option G_OPT_V_INPUT
#% key: calibration_points
-#% required: yes
+#%required: yes
#%end
#%option G_OPT_R_INPUT
#% key: additional_band1
-#% required: no
+#%required: no
#%end
#%option G_OPT_R_INPUT
#% key: additional_band2
-#% required: no
+#%required: no
#%end
#%option G_OPT_R_INPUT
#% key: additional_band3
-#% required: no
+#%required: no
#%end
#%option G_OPT_R_INPUT
#% key: additional_band4
-#% required: no
+#%required: no
#%end
#%option G_OPT_R_OUTPUT
#% key: depth_estimate
-#% required: yes
+#%required: yes
#%end
#%option
#% key: tide_height
-#% type: double
-#% multiple: no
-#% required: no
-#% description: Tide correction to the time of satellite image capture
+#%type: double
+#%multiple: no
+#%required: no
+#%description: Tide correction to the time of satellite image capture
#%end
#%option
#% key: calibration_column
-#% type: string
-#% multiple: no
-#% required: yes
-#% description: Name of the column which stores depth values
+#%type: string
+#%multiple: no
+#%required: yes
+#%description: Name of the column which stores depth values
#%end
#%flag
-#% key: b
+#%key: f
+#% description: select if only want to run Fixed-GWR model
+#%end
+#%flag
+#%key: b
#% description: select kernel function as bi-square
#%end
@@ -108,10 +111,18 @@
tide_height = options['tide_height']
calibration_column = options['calibration_column']
bisquare = flags['b']
- g.run_command('g.region', raster=Green)
+ fixed_GWR = flags['f']
+
+ res = g.parse_command('g.region', raster=Green, flags='g')
g.run_command('v.to.rast', input=Calibration_points, type='point',
use='attr', attribute_column=calibration_column,
output='tmp_Calibration_points')
+ # hull generation from calibration depth points
+ g.run_command('v.hull', input=Calibration_points, output='tmp_hull',
+ overwrite=True)
+ # buffer the hull to ceate a region for including all calibration points
+ g.run_command('v.buffer', input='tmp_hull', output='tmp_buffer',
+ distance=float(res['nsres']), overwrite=True)
if tide_height:
cal = g.parse_command('r.univar', map='tmp_Calibration_points',
flags='g')
@@ -136,17 +147,16 @@
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.message("Ditermining minimum value for %s" % i)
g.run_command('g.region', vector=Calibration_points)
- # Avoid zero values
+ # To ignore 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)
@@ -159,146 +169,183 @@
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.message("Deep water ditermination for %s" % i)
+ g.run_command('r.mask', vector='tmp_buffer', 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)
+ g.run_command('r.mask', raster='tmp_water', overwrite=True)
+ g.mapcalc("{tmp_crctd} = ({tmp_crctd} * 1)"
+ .format(tmp_crctd='tmp_crctd' + str(j)))
except:
g.message("Cannot find deep water pixels")
- g.run_command('g.remove', type='raster', name='MASK', flags='f')
+ g.run_command('r.mask', vector='tmp_buffer', overwrite=True)
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)
+ g.run_command('r.mask', raster='tmp_water', overwrite=True)
+ g.mapcalc("{tmp_crctd} = ({tmp_crctd} * 1)"
+ .format(tmp_crctd='tmp_crctd' + str(j)))
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-' \
+ try:
+ Module('r.gwr')
+ except:
+ g.run_command('g.extension', extension='r.gwr')
+ if fixed_GWR:
+ if bisquare:
+ g.message("Calculating optimal bandwidth using bisqare kernel...")
+ bw = g.parse_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points',
+ kernel='bisquare', flags='ge')
+ g.message("Running Fixed-GWR using bisqare kernel...")
+ g.run_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points',
+ estimates='tmp_bathymetry', kernel='bisquare',
+ bandwidth=int(bw['estimate']))
+ else:
+ g.message("Calculating optimal bandwidth using gaussian kernel...")
+ bw = g.parse_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points', flags='ge')
+ g.message("Running Fixed-GWR using gaussian kernel...")
+ g.run_command('r.gwr', mapx=crctd_lst,
+ mapy='tmp_Calibration_points',
+ estimates='tmp_bathymetry',
+ bandwidth=int(bw['estimate']))
+ else:
+ try:
+ # For GWmodel in R
+ 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,'
+ 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',
+ 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 ,' \
+ 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,' \
+ 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],' \
+ 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')
+ 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' \
+ 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" \
+ 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('\nCalculating optimal bandwidth using " \
"bisquare kernel..\n')\n"
- ref_flag += 'BW_Rapid_ref.sdf=bw.gwr(tmp_Calibration_points~%s,' \
+ 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' \
+ ref_flag += "cat('\nCalculating euclidean distance\n')\n"
+ 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' \
+ ref_flag += "cat('\nRunning A-GWR using bisquare kernel\n')\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,' \
+ r_file.write(ref_flag)
+ if not bisquare:
+ ref_fla = "cat('\nCalculating optimal bandwidth 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' \
+ ref_fla += "cat('\nCalculating euclidean distance\n')\n"
+ 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' \
+ ref_fla += "cat('\nRunning A-GWR using gaussian kernel\n')\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=",",' \
+ 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',
+ 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,
+ 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,
+ 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,
+ 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,
+ g.run_command('r.gwr', mapx=crctd_lst,
mapy='tmp_Calibration_points',
estimates='tmp_bathymetry',
bandwidth=int(bw['estimate']))
@@ -313,7 +360,7 @@
def cleanup():
g.run_command('g.remove', type='raster', pattern='MASK', flags='f')
- g.run_command('g.remove', type='raster', pattern='*tmp*', flags='f')
+ g.run_command('g.remove', type='all', pattern='*tmp*', flags='f')
try:
os.remove('crctd.R')
os.remove('prediction.txt')
@@ -323,4 +370,3 @@
if __name__ == '__main__':
atexit.register(cleanup)
main()
-
More information about the grass-commit
mailing list