[GRASS-dev] New module for lake filling - r.lake

Re-posting, as previous mail was larger than lists accept.

Hello all GRASSers.

GRASS CVS contains new module r.lake. This module can be used to create lake
maps if only place and lake level is known. To use this module You need
terrain raster map (DEM), water level for lake in same units as DEM values
and coordinates for any point in lake.

Module uses 3x3 sliding window to fill whole area connected to seed point (can
be coordinate pair or already existing raster map) and having height values
below water level. Output map will have lake depth in DEM units. More
information about module is in help page.

I use this module for palaeolake reconstructions based on single shore line
measurement. It also can be usefull for reservoir planing, creating maps with
scenarios for lake level changes etc. Unfortunately mailing lists does not
accept attachments larger than 100kb - no sample images with r.lake in
action.

Currently this module has some limitations and I ask for help to other GRASS
developers to help me to fix them.
Most important - module is NOT large file safe.
Lake volume calculations work only if DEM is in meters.
G_area_of_cell_at_row() always returns cell area in square meters, but I
could not figure out how to get information about DEM units. Thus if both are
in meters, volume is correct, else not ...

Thanks to all who helped me to code this module.

Any comments, ideas, enhancements - always welcome.

Maris Nartiss.

On Jun 17, 2006, at 5:26 AM, Māris Nartišs wrote:

Re-posting, as previous mail was larger than lists accept.

Hello all GRASSers.

GRASS CVS contains new module r.lake. This module can be used to create lake
maps if only place and lake level is known. To use this module You need
terrain raster map (DEM), water level for lake in same units as DEM values
and coordinates for any point in lake.

Module uses 3x3 sliding window to fill whole area connected to seed point (can
be coordinate pair or already existing raster map) and having height values
below water level. Output map will have lake depth in DEM units. More
information about module is in help page.

I use this module for palaeolake reconstructions based on single shore line
measurement. It also can be usefull for reservoir planing, creating maps with
scenarios for lake level changes etc. Unfortunately mailing lists does not
accept attachments larger than 100kb - no sample images with r.lake in
action.

Currently this module has some limitations and I ask for help to other GRASS
developers to help me to fix them.
Most important - module is NOT large file safe.
Lake volume calculations work only if DEM is in meters.
G_area_of_cell_at_row() always returns cell area in square meters, but I
could not figure out how to get information about DEM units. Thus if both are
in meters, volume is correct, else not ...

we have already discussed this in relation to r.slope.aspect which has the same problem,
currently it gives warning when x,y in feet is converted to meters so that user
can check and make sure that the elevation is in meters too.
There is no place currently in GRASS to store DEM z-units (and vertical datum) so it is left to
the user to make sure that the units are handled correctly (big pain here in US where
feet and meters are about equally common and mixed in various ways).

I believe that for GRASS7 we should include the vertical datum and units information for each location even if we
do not provide the transformation tools - this should help users better manage
their elevation and volume data and they can use external tools to convert the vertical datum
for the data that they want to import to the datum used in the GRASS location.
For example there is an official open source code for US coastal vertical datums conversion
available here (it is written in Java2 though)

http://chartmaker.ncd.noaa.gov/csdl/vdatum.htm

This would not prevent the users to mix misc. vert. and horiz. datums and units in their database
(at their own peril as I have learned the hard way), but it will encourage the management
of data in a well defined coordinate system and units in all 3 dimensions.

Helena

P.S. I would like to add this suggestion to the wiki development page - apparently it should go
under Plans but it does not fit into Overview, should I add it to sand box?
Or should we start a GRASS7 proposals/suggestions/todo section?

Thanks to all who helped me to code this module.

Any comments, ideas, enhancements - always welcome.

Maris Nartiss.

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

has anybody seen (or have an explanation) for the following vm_allocate error
that we got from v.proj when projecting rather large point file on Mac running recent CVS GRASS
(I was able to run v.proj with the same data on a FC5 machine without problem).
We got similar error from v.surf.rst after 86% done on 5000x5000 raster (but I was able
to compute bigger raster on a smaller mac laptop). I am trying to find out whether it is the data
or something specific to Mac (a 2GB memory machine) or a bug that is causing this error.

thanks, Helena

GRASS 6.1.cvs (jockeyspft):~ > v.proj input=earl030921oregoni
location=coastll mapset=PERMANENT output=earl030921oregoni
Input Projection Parameters: +proj=latlong +a=6378137 +rf=298.257222101
+no_defs +towgs84=0.000000,0.000000,0.000000
Input Unit Factor: 1
Output Projection Parameters: +proj=lcc +x_0=0.6096012192024384e+06 +y_0=0
+lon_0=79dw +lat_0=33d45'n +lat_1=36d10'n +lat_2=34d20'n +a=6378137
+rf=298.257222101 +no_defs +towgs84=0.000000,0.000000,0.000000
Output Unit Factor: 0.3048006096012192
Creating vector file...
Building topology ...
Registering lines: v.proj(2057) malloc: *** vm_allocate(size=8421376)
failed (error code=3)
v.proj(2057) malloc: *** error: can't allocate region
v.proj(2057) malloc: *** set a breakpoint in szone_error to debug
node.c:48: failed assertion `n'
Abort trap
--- end printout 030921-----

Although the vector was created, I wasn't able to display it. On the
contrary, I was able to display the 030916 vector despite reporting a
similar error without abortion:
---- error printout 030916---
(v.proj(753) malloc: *** error: can't allocate region
v.proj(753) malloc: *** set a breakpoint in szone_error to debug
v.proj(753) malloc: *** vm_allocate(size=24088576) failed (error code=3)
---- end error printout 030916---

there may be a huge segment in the masked area, but the error message is not comming from the rst code
v.surf.rst input=L18 layer== dmax=3.0 dmin=1.0 elev=surf500501Lidar maskmap=MASK zmult=1.0 tension=40.0
smooth=10.0 segmax=20 npmin=120

   86%
WARNING: taking too long to find points for interpolation--please change
          the region to area where your points are. Continuing
          calculations...
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
ERROR: G_malloc: out of memory

Helena Mitasova
Dept. of Marine, Earth and Atm. Sciences
1125 Jordan Hall, NCSU Box 8208,
Raleigh NC 27695
http://skagit.meas.ncsu.edu/~helena/

On Jun 28, 2006, at 8:24 PM, Helena Mitasova wrote:

On Jun 17, 2006, at 5:26 AM, Māris Nartišs wrote:

Re-posting, as previous mail was larger than lists accept.

Hello all GRASSers.

GRASS CVS contains new module r.lake. This module can be used to create lake
maps if only place and lake level is known. To use this module You need
terrain raster map (DEM), water level for lake in same units as DEM values
and coordinates for any point in lake.

Module uses 3x3 sliding window to fill whole area connected to seed point (can
be coordinate pair or already existing raster map) and having height values
below water level. Output map will have lake depth in DEM units. More
information about module is in help page.

I use this module for palaeolake reconstructions based on single shore line
measurement. It also can be usefull for reservoir planing, creating maps with
scenarios for lake level changes etc. Unfortunately mailing lists does not
accept attachments larger than 100kb - no sample images with r.lake in
action.

Currently this module has some limitations and I ask for help to other GRASS
developers to help me to fix them.
Most important - module is NOT large file safe.
Lake volume calculations work only if DEM is in meters.
G_area_of_cell_at_row() always returns cell area in square meters, but I
could not figure out how to get information about DEM units. Thus if both are
in meters, volume is correct, else not ...

we have already discussed this in relation to r.slope.aspect which has the same problem,
currently it gives warning when x,y in feet is converted to meters so that user
can check and make sure that the elevation is in meters too.
There is no place currently in GRASS to store DEM z-units (and vertical datum) so it is left to
the user to make sure that the units are handled correctly (big pain here in US where
feet and meters are about equally common and mixed in various ways).

I believe that for GRASS7 we should include the vertical datum and units information for each location even if we
do not provide the transformation tools - this should help users better manage
their elevation and volume data and they can use external tools to convert the vertical datum
for the data that they want to import to the datum used in the GRASS location.
For example there is an official open source code for US coastal vertical datums conversion
available here (it is written in Java2 though)

http://chartmaker.ncd.noaa.gov/csdl/vdatum.htm

This would not prevent the users to mix misc. vert. and horiz. datums and units in their database
(at their own peril as I have learned the hard way), but it will encourage the management
of data in a well defined coordinate system and units in all 3 dimensions.

Helena

P.S. I would like to add this suggestion to the wiki development page - apparently it should go
under Plans but it does not fit into Overview, should I add it to sand box?
Or should we start a GRASS7 proposals/suggestions/todo section?

Thanks to all who helped me to code this module.

Any comments, ideas, enhancements - always welcome.

Maris Nartiss.

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

Brad Douglas wrote:

Helena,

Those are error messages from the operating system kernel. It has
literally run out of RAM (or it thinks it has - e.g. 32bit coding on
64bit system).

I suppose it may depend on where it runs out of memory in the file as to
whether the results are displayable or not (but I would expect it to be
incomplete, regardless).

I didn't realize v.proj was such a demanding memory hog.
  

Brad, thank you, this helps. When you look at the output it is not v.proj itself that is the problem but
the usual topology building - so it converts the points and then builds the topology and that is when
it runs out of memory - a well known problem. Except that this is a little bit nicer behavior than the usual
machine freeze that you get when building topology on a large file.
With v.surf.rst it looks like we have a problem because there is no reason it should be running out of memory
when computing the segments, hopefully I will be able to find and fix that one but the topology stuff
seems to be beyond what anybody can do for now.

Helena

On Wed, 2006-06-28 at 22:38 -0400, Helena Mitasova wrote:
  

has anybody seen (or have an explanation) for the following vm_allocate error
that we got from v.proj when projecting rather large point file on Mac running recent CVS GRASS
(I was able to run v.proj with the same data on a FC5 machine without problem).
We got similar error from v.surf.rst after 86% done on 5000x5000 raster (but I was able
to compute bigger raster on a smaller mac laptop). I am trying to find out whether it is the data
or something specific to Mac (a 2GB memory machine) or a bug that is causing this error.

thanks, Helena

GRASS 6.1.cvs (jockeyspft):~ > v.proj input=earl030921oregoni
location=coastll mapset=PERMANENT output=earl030921oregoni
Input Projection Parameters: +proj=latlong +a=6378137 +rf=298.257222101
+no_defs +towgs84=0.000000,0.000000,0.000000
Input Unit Factor: 1
Output Projection Parameters: +proj=lcc +x_0=0.6096012192024384e+06 +y_0=0
+lon_0=79dw +lat_0=33d45'n +lat_1=36d10'n +lat_2=34d20'n +a=6378137
+rf=298.257222101 +no_defs +towgs84=0.000000,0.000000,0.000000
Output Unit Factor: 0.3048006096012192
Creating vector file...
Building topology ...
Registering lines: v.proj(2057) malloc: *** vm_allocate(size=8421376)
failed (error code=3)
v.proj(2057) malloc: *** error: can't allocate region
v.proj(2057) malloc: *** set a breakpoint in szone_error to debug
node.c:48: failed assertion `n'
Abort trap
--- end printout 030921-----

Although the vector was created, I wasn't able to display it. On the
contrary, I was able to display the 030916 vector despite reporting a
similar error without abortion:
---- error printout 030916---
(v.proj(753) malloc: *** error: can't allocate region
v.proj(753) malloc: *** set a breakpoint in szone_error to debug
v.proj(753) malloc: *** vm_allocate(size=24088576) failed (error code=3)
---- end error printout 030916---

there may be a huge segment in the masked area, but the error message is not comming from the rst code
v.surf.rst input=L18 layer== dmax=3.0 dmin=1.0 elev=surf500501Lidar maskmap=MASK zmult=1.0 tension=40.0
smooth=10.0 segmax=20 npmin=120

   86%
WARNING: taking too long to find points for interpolation--please change
          the region to area where your points are. Continuing
          calculations...
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
v.surf.rst(1340) malloc: *** vm_allocate(size=8421376) failed (error code=3)
v.surf.rst(1340) malloc: *** error: can't allocate region
v.surf.rst(1340) malloc: *** set a breakpoint in szone_error to debug
ERROR: G_malloc: out of memory

Helena Mitasova
Dept. of Marine, Earth and Atm. Sciences
1125 Jordan Hall, NCSU Box 8208,
Raleigh NC 27695
http://skagit.meas.ncsu.edu/~helena/

On Jun 28, 2006, at 8:24 PM, Helena Mitasova wrote:

On Jun 17, 2006, at 5:26 AM, Māris Nartišs wrote:

Re-posting, as previous mail was larger than lists accept.

Hello all GRASSers.

GRASS CVS contains new module r.lake. This module can be used to create lake
maps if only place and lake level is known. To use this module You need
terrain raster map (DEM), water level for lake in same units as DEM values
and coordinates for any point in lake.

Module uses 3x3 sliding window to fill whole area connected to seed point (can
be coordinate pair or already existing raster map) and having height values
below water level. Output map will have lake depth in DEM units. More
information about module is in help page.

I use this module for palaeolake reconstructions based on single shore line
measurement. It also can be usefull for reservoir planing, creating maps with
scenarios for lake level changes etc. Unfortunately mailing lists does not
accept attachments larger than 100kb - no sample images with r.lake in
action.

Currently this module has some limitations and I ask for help to other GRASS
developers to help me to fix them.
Most important - module is NOT large file safe.
Lake volume calculations work only if DEM is in meters.
G_area_of_cell_at_row() always returns cell area in square meters, but I
could not figure out how to get information about DEM units. Thus if both are
in meters, volume is correct, else not ...
        

we have already discussed this in relation to r.slope.aspect which has the same problem,
currently it gives warning when x,y in feet is converted to meters so that user
can check and make sure that the elevation is in meters too.
There is no place currently in GRASS to store DEM z-units (and vertical datum) so it is left to
the user to make sure that the units are handled correctly (big pain here in US where
feet and meters are about equally common and mixed in various ways).

I believe that for GRASS7 we should include the vertical datum and units information for each location even if we
do not provide the transformation tools - this should help users better manage
their elevation and volume data and they can use external tools to convert the vertical datum
for the data that they want to import to the datum used in the GRASS location.
For example there is an official open source code for US coastal vertical datums conversion
available here (it is written in Java2 though)

http://chartmaker.ncd.noaa.gov/csdl/vdatum.htm

This would not prevent the users to mix misc. vert. and horiz. datums and units in their database
(at their own peril as I have learned the hard way), but it will encourage the management
of data in a well defined coordinate system and units in all 3 dimensions.

Helena

P.S. I would like to add this suggestion to the wiki development page - apparently it should go
under Plans but it does not fit into Overview, should I add it to sand box?
Or should we start a GRASS7 proposals/suggestions/todo section?

Thanks to all who helped me to code this module.

Any comments, ideas, enhancements - always welcome.

Maris Nartiss.
        
--
Helena Mitasova
Department of Marine, Earth and Atmospheric Sciences
North Carolina State University
1125 Jordan Hall
NCSU Box 8208
Raleigh, NC 27695-8208
http://skagit.meas.ncsu.edu/~helena/

email: hmitaso@unity.ncsu.edu
ph: 919-513-1327 (no voicemail)
fax 919 515-7802

On Wed, 28 Jun 2006 22:38:12 -0400
Helena Mitasova <hmitaso@unity.ncsu.edu> wrote:

has anybody seen (or have an explanation) for the following
vm_allocate error
that we got from v.proj when projecting rather large point file on
Mac running recent CVS GRASS
(I was able to run v.proj with the same data on a FC5 machine
without problem).
We got similar error from v.surf.rst after 86% done on 5000x5000
raster (but I was able
to compute bigger raster on a smaller mac laptop). I am trying to
find out whether it is the data
or something specific to Mac (a 2GB memory machine) or a bug that is
causing this error.

workaround:

use m.proj (or just cs2cs, which m.proj is a front-end for) to reproject
the data and then v.in.ascii to import the data into the target
location.

you can even pipe things:

# assume LL WGS84 input, "tr" and fs= as needed
cat xyz.txt | m.proj -i | v.in.ascii -bt out=xyz_proj

no idea if there is any memory bottlenecks using that method.

Hamish