[GRASS-user] Help needed with interactive selection of vectors

[Hamish]

Thomas Adams wrote:
> > I work for the US National Weather Service at the Ohio River 
> > Forecast Center. We generate, on a daily basis, a graphic that 
> > looks like this:
> > http://www.erh.noaa.gov/er/ohrfc/fop.html

the site seems unreachable today for me.. endless wait.

> > The linked graphic would be pretty trivial to create in GRASS using
> > ps.map. However, during times where there is a reasonable risk for
> > significant flooding, the map must be altered to designate the 
> > smaller subbasins (or areas) that are at risk and to annotate the 
> > period over
> > which the threat exists. I'm sorry I can not link an example to 
> > show, since there is no current flooding in the US or risk thereof.

it is simple to script ps.map with variable text & queries.
see example script below.

> > Currently, the process is carried out interactively using ArcGIS. 
> > I am trying to offer up an open source alternative with GRASS. I
> > am sure a
> > more sensible, non-interactive way of generating the map explicitly
> > through a modeling approach would be best (and we do, of course, 
> > base designating the basins at risk on hydrologic modeling), but 
> > there is no *direct* link. Currently human interaction is used and
> > some will insist on this.

Human interaction is fine, at times of public risk human oversight is critical.
Automation is a great help though, and if the digitizing work can be 95%
prepared for the operator - and they just have to touch it up - then all the
better for speed of deployment & reduction of opportunities for human error.

> > So, I am left wondering how to best approach this is GRASS, 
> > namely, either:
> >
> > (1) have the ability to manually select smaller subbasins and 
> > group them together to designate the area at risk, or

see script below for generating map simply from specified catchment IDs.

> > (2) hand draw/digitize the area

Digitize areas with v.digit over backdrop image. Backdrop image can be
static or with raster overlay from hydrologic model output.

(or digitize with QGIS, or digitizing frontend like OpenLayers WPS- try the
polygon tool: http://www.les-ejk.cz/english/pywps-and-openlayers )


[...]
> Sorry for not being clear. Basically now the mode of operation is:
> 
> (1) launch ArcGIS which has been customized somewhat to do the following
> (some background coverages appear automatically for reference, as the link
> shows)...
> (2) the user manually digitizes polygonal areas corresponding to the areas of
> flood risk

this seems highly inefficient to me.

> (3) the polygons are annotated to show the period over which they are valid
> (4) the custom control buttons allow the user to export the map and to have a
> text file created and have it ftp'ed to another location
> 
> Now, all of this is pretty much doable in GRASS, though the GUI is not as
> refined and some coding changes would be needed (I think) to 'duplicate' this
> process.

It seems best done from a script, in which case the GUI is moot and the CLI is
king. If partial interaction is required, the script could be written to pop up
the digitizing tool as needed, then preview, "Accept or Re-edit", and finally
continue the script to do the rest.

A g.parser script could be made friendly & have checkboxes (e.g. the v.clean
GUI tool selection) to activate catchment by names/ID. Or a fancy front end
could allow the user to click on an area of a map. (e.g. parse 'd.what.vect
-xft'; or generate from the HTMLMAP monitor driver [see d.mon]; or pyWPS style
solution)

> I'm trying to avoid this right now for several reasons, including
> time constraints. What I thought *could* be done, would be to have some
> reference subbasin boundaries displayed in the GRASS display window along
> with the other reference vector maps (state political boundaries and our area
> of responsibility). Allow the user to interactively select the basins at risk
> and have those areas ONLY display as areas of risk.

hopefully the model output will tell you (ie the script) which subbasins are at
risk...


> Then the process would be similar to (1)—(4), except that:
> 
> (1) I do not see how polygons can be manually selected and have the option to
> remove unselected basins/polygons and

limit by known catchment ID (vector feature category) or SQL query on vector
area attribute (catchment's name). Make the query part of the script set by a
environment variable, see example below.

> (2) I do not see how to manually annotate a polygon in such a way that the
> polygon(s) are clearly tied to the annotation.

In ps.map there is the vlegend instruction which does this automatically.
d.legend will display raster category labels, if present, and hide cat numbers.

> However, this may all be moot because, it seems, a significant number of
> people want to have the ability to freehand draw the areas of risk, so that
> the selection process I was suggesting is not of interest.

Again, preparing a proto-map for the human digitizer to finish/check may be of
great value, both by speeding their workflow and reducing chances for error.
Manual-override at key steps would be very important of course.


Here is a little demo ps.map script I've thrown together demonstrating the
idea.
It builds on the Spearfish catchments generated in the (brand new) r.watershed
help page example*. All you need to specify is catchments at risk, it does the
rest.

[*]
http://freegis.org/cgi-bin/viewcvs.cgi/*checkout*/grass6/raster/r.watershed/front/description.html?rev=HEAD&content-type=text/html

screenshots of output for whole catchment warning and lowlying areas warning:
  http://bambi.otago.ac.nz/hamish/grass/screenshots/flood_warning_basins.png
  http://bambi.otago.ac.nz/hamish/grass/screenshots/flood_warning_lowlying.png
screenshot of Spearfish catchments:
  http://grass.itc.it/screenshots/raster.php


First we need to generate polygons for low lying ares from each catchment
created with the r.watershed help page example:

# create danger polygons for each catchment
#  (homebrew example of predictor for at-risk lowlying areas)
r.cost -k in=slope out=slope_cost_from_rivers start_rast=rwater.course
r.colors slope_cost_from_rivers col=bcyr -e

# choose some threshold. this could be based on hydrologic model output.
#  for now just set for lower half of slope-costs in the map
eval `r.univar -g slope_cost_from_rivers`
r.mapcalc "flood.zones = if(slope_cost_from_rivers < $mean, \
     rwater.basin, null() )"

# convert to vector areas
r.to.vect -s in=flood.zones out=flood_zones feature=area
v.db.dropcol map=flood_zones column=label
v.db.renamecol map=flood_zones column=value,catchment



################################
# dynamic bits, set by model or interactive selection
#  (eg scripted 'd.what.* -xft' to highlight areas)
# Here we say catchment #6 and #12 are at risk.

FLOOD_BASINS=6,12
TIMESTAMP=`date`



###############################
# the rest is a static template
#  commands could be written to a file and then passed to ps.map input=

# convert catchment list to SQL query
FLOOD_BASIN_SQL=`echo $FLOOD_BASINS | \
  sed -e 's/^/catchment = /' -e 's/,/ OR catchment = /g'`

# generate PostScript image
ps.map out=flood_warning.ps << EOF
raster elevation.dem
# example of dashed lines
vlines roads
  style 000011
  where label ~ 'highway'
  width 0.25
  label Highways
  end
vlines rwater_course
  color grey
  label Water courses
  end
vareas rwater_basins
#change to "flood_zones" vector map for lowlying areas instead of
#  entire catchment
  where $FLOOD_BASIN_SQL
  color red
  fcolor red
  width 2
  pat $GISBASE/etc/paint/patterns/diag_up.eps
  label Catchements subject to flood risk
  end
vlegend
  end
text 100% -3% Warning issued $TIMESTAMP
  font Helvetica-BoldOblique
  fontsize 12
  ref right
  end
end
EOF


# convert to PDF
ps2pdf flood_warning.ps
# convert to PNG
pstoimg -antialias -aaliastext -out flood_warning.png -scale 1.3 \
  -crop tb flood_warning.ps


The PDF output looks the best.


Rendering problems:
* EPS hatch pattern renders in gv and PNG with black boxes behind lines??
    These do not appear in PDF + acrobat reader.
* EPS hatch patterns are tiled so fill lines have slightly noticeable breaks.
* We need to add something to allow skipping vectors in ps.map's vlegends
    ideas: labelskip [y|N],  test for label='' or 'none' or 'skip'?



Hamish


__________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around 
http://mail.yahoo.com