[GRASS-user] r.fillnulls for large area (Nick Jachowski)

Dear Nick,

I am working on SRTM rasters of resolution 1km, 500m, 250m and recently 90m using GRASS for long and had countered the problem being faced by you.

My workaround is given below and as a by-product you may get a very nice coastline vector/raster also.

  1. I have converted all the nulls with r.mapcalc to -160 using the formula C=if(A, if(isnull(B),-160,B))
    Where,
    A: region of interest
    B: original srtm raster with nulls for sea and holes
    C: Output map in which nulls have been converted to -160.

(Why -160, because in my map, as from my human memory, the lowest value was -70. I just chose a value lower than this value. We will need it later on to convert back them as holes for r.fillnulls to work upon.)

A sample from the history file of the raster is produced below:

Sun Aug 1 17:49:10 2010
test5
PERMANENT
pks
raster

generated by r.mapcalc
if(test_nL_region@PERMANENT, if(isnull(ib2@PERMANENT), -160, ib2@PERMANENT))

A: test_nL_region

This raster map is a region of my interest out the big srtm raster converted from a vector. Grass has a command , which converts the current region to a vector file.
Again, the sample from the history file of the raster is produced below:

ed Jul 28 10:09:42 2010
test_nL_region
PERMANENT
pks
raster
Vector Map: nL_region@PERMANENT in mapset PERMANENT
Original scale from vector map: 1:1
generated by v.to.rast
v.to.rast input=“nL_region@PERMANENT” output=“test_nL_region” use=“c
at” type=“point,line,area” layer=1 value=1 rows=4096

nL_region:
Again, the sample from the history file of the vector is produced below:
COMMAND: v.in.region output=“nL_region” type=“area” cat=1
GISDBASE: /home/pks/renamed_grassdata
LOCATION: newLocation MAPSET: PERMANENT USER: pks DATE: Wed Jul 28 09:53:43 2010

So, now we have -160s all around , i.e. the sea and holes on land.

Now use the r.reclass.area feature to differentiate sea from holes on land. A sample from my recent work is given below:

Tue Nov 8 13:55:47 2011
bi_reclassed3_area.recl
work2
rdcgis1
reclass
Reclassified map based on:
Map [bi_reclassed3.clump.bi_reclassed3_area] in mapset [work2]
generated by r.reclass

I am using GRASS7 and in the history given above , all the parameters are not available.
In my case , I used 1 billion hectare as “greater than area” parameter to bring out the sea , land , and holes in land in different categories.

Thereafter, use mask and pick up the landmass with holes in land , which you want to fill.
Just remember to convert back the -160 to null again.

I felt that r.fillnulls works best when the nulls have values around them.

By the way, CGIAR provides holes-filled-SRTM-raster.You may like to read their literature for understanding the algorithms used to fill the nulls. It’s good and useful.

Try this link for downloading clean rasters of your area of interest:

http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

Hope, it helps.

Some of my observations are listed below:

  1. If the system lacks RAM, then increase patience.
  2. If the system lacks hard-disk space, don’t go for GIS.
  3. Learn python.

On Fri, Dec 2, 2011 at 5:37 AM, <grass-user-request@lists.osgeo.org> wrote:

Send grass-user mailing list submissions to
grass-user@lists.osgeo.org

To subscribe or unsubscribe via the World Wide Web, visit
http://lists.osgeo.org/mailman/listinfo/grass-user
or, via email, send a message with subject or body ‘help’ to
grass-user-request@lists.osgeo.org

You can reach the person managing the list at
grass-user-owner@lists.osgeo.org

When replying, please edit your Subject line so it is more specific
than “Re: Contents of grass-user digest…”

Today’s Topics:

  1. Re: r.fillnulls for large area (Nick Jachowski)
  2. topology model resume (and some proposals) (G. Allegri)
  3. Re: [GRASS-dev] Should v.in.ogr clean topology by default ?
    (Hamish)
  4. Re: r.fillnulls for large area (Hamish)
  5. Re: r.fillnulls for large area (Marcello Gorini)

Message: 1
Date: Fri, 2 Dec 2011 00:07:07 +0700
From: Nick Jachowski <njachowski@gmail.com>
Subject: Re: [GRASS-user] r.fillnulls for large area
To: Hamish <hamish_b@yahoo.com>
Cc: grass-user@lists.osgeo.org
Message-ID:
<CAPFL4TFhef32rraoch_+29xet8xWm6-B4sWYtva1qKi4GDC7qQ@mail.gmail.com>
Content-Type: text/plain; charset=“iso-8859-1”

Hi Hamish,

Thanks for the advice to fill the ocean with zeros. It worked better, and
didn’t get “killed”, but it still didn’t work…This time it made it
further than before, but it got “abondoned” - see the output below…

Maybe I need to split my region into subregions and run r.fillnulls on each
subregion then patch them together. If anybody has a better idea I would be
grateful to know.

Thanks!
Nick

GRASS 6.4.1 (seasia_ll):~ > r.fillnulls in=srtmm0 out=srtmm0f
Locating and isolating NULL areas…
100%
Reading input raster map r_fillnulls_18931@PERMANENT
100%
Finding buffer zones…
100%
Writing output raster map <r_fillnulls_18931.buf>…
100%
100%
Creating interpolation points…
Extracting points…
100%
Building topology for vector map <vecttmp_fillnulls_18931>…
Registering primitives…
13428068 primitives registered
13428068 vertices registered
Building areas…
100%
0 areas built
0 isles built
Attaching islands…
Attaching centroids…
100%
Number of nodes: 13428068
Number of primitives: 13428068
Number of points: 13428068
Number of lines: 0
Number of boundaries: 0
Number of centroids: 0
Number of areas: 0
Number of isles: 0
ERROR: /usr/lib/grass64/scripts/r.fillnulls abandoned. Removing temporary
maps, restoring user mask if needed:
Removing raster
Removing raster <r_fillnulls_18931>
Removing raster <r_fillnulls_18931.buf>
Removing raster <r_fillnulls_18931_filled>
WARNING: Raster map <r_fillnulls_18931_filled> not found
WARNING: <r_fillnulls_18931_filled> nothing removed
Removing vector <vecttmp_fillnulls_18931>
WARNING: Table <vecttmp_fillnulls_18931> linked to vector map
<vecttmp_fillnulls_18931> does not exist
WARNING: raster <usermask_mask.18931> not found

On Thu, Dec 1, 2011 at 8:36 AM, Hamish <hamish_b@yahoo.com> wrote:

Nick wrote:

I am having trouble filling the null areas in the SRTM (3 second
resolution) DEM for Southeast Asia. I think the reason is because the
region I am trying fill is too large. The region is almost a billion
cells - see g.region -p output below - but about 40% of that region is
ocean areas with no DEM data. I’m wondering if maybe the problem is that
r.fillnulls is trying to fill the ocean areas.

it will try that. a good first step is to fill in the ocean with 0s.
I was just doing this yesterday with srtm data actually.
GRASS> r.mapcalc “zero = 0”

then r.patch with your srtm data listed first and the zeros listed last.

but be careful you don’t fill in zeros where there are holes in the data
on land. if you have a vector coastline use “v.to.rast use=val val=1” to
create a land mask, then invert with:
r.mapcalc “sea.mask = if(isnull(land.mask), 1, null())”

and use the sea.mask as the background fill-in for the r.patch step.

r.fillnulls gets to the point where it says “Building areas…” then
it says “Killed” and prints an error message.

often “Killed” happens when the computer ran out of memory.
not interpolating over the ocean may make the thing go much faster and
need much less memory, as it has much less to do.

I am using GRASS 6.4.1. Eventually I plan on using the null-filled
DEM to delineate the watersheds in Southeast Asia.

start thinking about your watershed threshold size now :slight_smile:

good luck,
Hamish

-------------- next part --------------
An HTML attachment was scrubbed…
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20111202/2a3d50ef/attachment-0001.html


Message: 2
Date: Thu, 1 Dec 2011 18:33:01 +0100
From: “G. Allegri” <giohappy@gmail.com>
Subject: [GRASS-user] topology model resume (and some proposals)
To: grass-user@lists.osgeo.org
Message-ID:
<CAB4g1=xqBcQOjoZ49n0adS0HGsEokccwtpS4nNtotzfMgKLvSQ@mail.gmail.com>
Content-Type: text/plain; charset=“iso-8859-1”

I resume (first as a repeat to myself) what I’ve learned from the various
email on the topic

Vectors can be:

LEVEL 1:

  • no topology → very limited use
    LEVEL 2:
  • unclean topology → limited use
  • clean topology → full support

I previously thought that LEVEL 2 was only possible for clean topologies,
and I was wrong…

At the moment there isn’t a tool to list the the uncorrect geometries from
a topological point of view. v.build only checks some constraints, not all.
The proposal is to extend it to check against all the rules that are
required to consider a geometry topologically correct (an extended flag to
v defaul.build maybe).

v.in.ogr builds and cleans (by default). It would be useful to have the
“clean” phase available to be launched independently. I mean, something
like an “automatic” flag for v.clean, that would operate the same cleaning
as during the import of a vector.

Conclusions: the topological correctness isn’t a constraint for the vector
topology data structure. GRASS haven’t all the topology rules hard-coded
(… or yes?). Most of thems (all?) are defined inside the code of v.build
and v.clean, but I suppose that there isn’t an autonomous
library/functionality that provide the semantics of a “correct topology”.
Am I wrong?

Thanks everyone for the support :wink:
giovanni
-------------- next part --------------
An HTML attachment was scrubbed…
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20111201/48197e88/attachment-0001.html


Message: 3
Date: Thu, 1 Dec 2011 13:01:38 -0800 (PST)
From: Hamish <hamish_b@yahoo.com>
Subject: [GRASS-user] Re: [GRASS-dev] Should v.in.ogr clean topology
by default ?
To: grass-user@lists.osgeo.org
Message-ID:
<1322773298.49418.YahooMailClassic@web110005.mail.gq1.yahoo.com>
Content-Type: text/plain; charset=us-ascii

[cc’d from grass-dev ML]

Hi,

a couple of points–

  • Radim said many times: v.in.ogr’s cleaning and v.clean are not
    the same thing. (with a number of !! added) He explained it and
    understood it always better than I ever could, so check in the
    mailing list archives for details.

  • GRASS is a topological GIS. so when data is imported into GRASS
    it is natural and desirable that it gets converted into the
    family way at that time.

  • In defense of the option to allow non topological imports:
    in the past I have imported a shapefile created by leading
    proprietary non-topological GIS. this file was a string of
    buffers around points with unique ID numbers. the idea was to
    run (the predecessor to) v.rast.stats on each buffer to collect
    some neighborhood context stats for each ID point, which could
    be loaded into a table. wherever the line of points was not
    exactly straight the buffers overlapped a little and grass’s
    cleaning would treat it as a ven diagram and create a new
    centroid in the overlapping sliver. no good for extracting the
    buffer area where=‘NODE_ID = 1234’.
    also you might consider a vector area map with multiple over-
    lapping home ranges for different territorial animals. not 3d,
    but still merging/flattening overlapping areas is not appropriate.

for things like neighboring land parcels, political lines, or
land/sea boundaries non-topological is a nightmare, but for
other tasks it can occasionally be very useful.

Hamish


Message: 4
Date: Thu, 1 Dec 2011 13:16:45 -0800 (PST)
From: Hamish <hamish_b@yahoo.com>
Subject: Re: [GRASS-user] r.fillnulls for large area
To: Nick Jachowski <njachowski@gmail.com>
Cc: grass-user@lists.osgeo.org
Message-ID:
<1322774205.78253.YahooMailClassic@web110013.mail.gq1.yahoo.com>
Content-Type: text/plain; charset=iso-8859-1

Nick wrote:

Thanks for the advice to fill the ocean with
zeros. It worked better, and didn’t get “killed”,
but it still didn’t work…This time it made it
further than before, but it got “abondoned” - see
the output below…

try editing the r.fillnulls script, on line 165
replace “r.to.vect input=…” with
“r.to.vect -b input=…”. I think it would work,
but I’m not sure if it would make it slower to run
v.surf.rst or not. (?)

Maybe I need to split my region into subregions
and run r.fillnulls on each subregion then patch
them together.

yes, that should work well.

Number of nodes: 13428068

building topology for anything more that 3-5 million
cells will consume a large amount of memory due to
the cumulative effect of the small, but real,
memory requirement of each individual feature.
maybe if you had 16gb RAM 13 million would be ok?

Hamish


Message: 5
Date: Thu, 1 Dec 2011 20:04:38 -0200
From: Marcello Gorini <gorini@gmail.com>
Subject: Re: [GRASS-user] r.fillnulls for large area
To: Hamish <hamish_b@yahoo.com>
Cc: grass-user@lists.osgeo.org
Message-ID:
<CAEvwxEJtmTd_dEfowfDcpfk7trjAWC=OzC7xjMZKQyPO4zp3Fg@mail.gmail.com>
Content-Type: text/plain; charset=“iso-8859-1”

Nick wrote:

Maybe I need to split my region into subregions
and run r.fillnulls on each subregion then patch
them together.

Hamish:

yes, that should work well.

I always wanted to create a script to do that. Since (I think) r.fill.nulls
identifies the null locations and creates small buffers around them, I
believe that if we combined something like the region handling of, say, the
r.roughness script, with r.fill.nulls, that would really speed up matters.

Don’t you think?

Cheers,
Marcello.
-------------- next part --------------
An HTML attachment was scrubbed…
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20111201/02e8ea3f/attachment.html



grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

End of grass-user Digest, Vol 68, Issue 2