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

Pankaj Kr Sharma pkscwc at gmail.com
Sat Dec 3 02:41:07 EST 2011


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 at PERMANENT, if(isnull(ib2 at PERMANENT), -160, ib2 at 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 at PERMANENT in mapset PERMANENT
Original scale from vector map: 1:1
generated by v.to.rast
v.to.rast input="nL_region at 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 at lists.osgeo.org> wrote:

> Send grass-user mailing list submissions to
>        grass-user at 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 at lists.osgeo.org
>
> You can reach the person managing the list at
>        grass-user-owner at 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 at gmail.com>
> Subject: Re: [GRASS-user] r.fillnulls for large area
> To: Hamish <hamish_b at yahoo.com>
> Cc: grass-user at lists.osgeo.org
> Message-ID:
>        <CAPFL4TFhef32rraoch_+29xet8xWm6-B4sWYtva1qKi4GDC7qQ at 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 at 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 <MASK>
> 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 at 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 :)
> >
> >
> > 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 at gmail.com>
> Subject: [GRASS-user] topology model resume (and some proposals)
> To: grass-user at lists.osgeo.org
> Message-ID:
>        <CAB4g1=xqBcQOjoZ49n0adS0HGsEokccwtpS4nNtotzfMgKLvSQ at 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 ;)
> 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 at yahoo.com>
> Subject: [GRASS-user] Re: [GRASS-dev] Should v.in.ogr clean topology
>        by      default ?
> To: grass-user at lists.osgeo.org
> Message-ID:
>        <1322773298.49418.YahooMailClassic at 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 at yahoo.com>
> Subject: Re: [GRASS-user] r.fillnulls for large area
> To: Nick Jachowski <njachowski at gmail.com>
> Cc: grass-user at lists.osgeo.org
> Message-ID:
>        <1322774205.78253.YahooMailClassic at 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 at gmail.com>
> Subject: Re: [GRASS-user] r.fillnulls for large area
> To: Hamish <hamish_b at yahoo.com>
> Cc: grass-user at lists.osgeo.org
> Message-ID:
>        <CAEvwxEJtmTd_dEfowfDcpfk7trjAWC=OzC7xjMZKQyPO4zp3Fg at 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 at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
>
>
> End of grass-user Digest, Vol 68, Issue 2
> *****************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20111203/76270677/attachment-0001.html


More information about the grass-user mailing list