[GRASS-user] r.fillnulls for large area

Nick_Jachowski · November 30, 2011, 4:07pm

Hello GRASS community,

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.

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

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

Anybody know what I’m doing wrong?

Thanks,
Nick

output from r.fillnulls:

Locating and isolating NULL areas…
100%
Reading input raster map r_fillnulls_11665@PERMANENT…
100%
Finding buffer zones…
100%
Writing output raster map <r_fillnulls_11665.buf>…
100%
100%
Creating interpolation points…
Extracting points…
100%
Building topology for vector map <vecttmp_fillnulls_11665>…
Registering primitives…
14602532 primitives registered
14602532 vertices registered
Building areas…
Killed
ERROR: /usr/lib/grass64/scripts/r.fillnulls abandoned. Removing temporary
maps, restoring user mask if needed:
Removing raster
Removing raster <r_fillnulls_11665>
Removing raster <r_fillnulls_11665.buf>
Removing raster <r_fillnulls_11665_filled>

GRASS 6.4.1 (seasia_ll):/media/little/projects/seasia > g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 31:00:01.5N
south: 0:58:31.5N
west: 88:59:58.5E
east: 112:00:01.5E
nsres: 0:00:03
ewres: 0:00:03
rows: 36030
cols: 27601
cells: 994464030

M_S · November 30, 2011, 6:56pm

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.

Maybe apply a raster mask over the ocean area?

Mark

Nick_Jachowski · November 30, 2011, 10:51pm

Hi Mark,

I tried that, and it made it farther (almost to completion I think), but it still died. See output below.

Thanks,
Nick

GRASS 6.4.1 (seasia_ll):/media/little/projects/seasia > r.fillnulls srtmm out=srtmmf
A user raster mask (MASK) is present. Saving it…
Rename raster to <usermask_mask.14798>
Locating and isolating NULL areas…
100%
Reading input raster map r_fillnulls_14798@PERMANENT…
100%
Finding buffer zones…
100%
Writing output raster map <r_fillnulls_14798.buf>…
100%
100%
Creating interpolation points…
Extracting points…
100%
Building topology for vector map <vecttmp_fillnulls_14798>…
Registering primitives…
13684301 primitives registered
13684301 vertices registered
Building areas…
100%
0 areas built
0 isles built
Attaching islands…
Attaching centroids…
100%
Killed
ERROR: /usr/lib/grass64/scripts/r.fillnulls abandoned. Removing temporary
maps, restoring user mask if needed:
Removing raster
Removing raster <r_fillnulls_14798>
Removing raster <r_fillnulls_14798.buf>
Removing raster <r_fillnulls_14798_filled>
WARNING: Raster map <r_fillnulls_14798_filled> not found
WARNING: <r_fillnulls_14798_filled> nothing removed
Removing vector <vecttmp_fillnulls_14798>
WARNING: Table <vecttmp_fillnulls_14798> linked to vector map
<vecttmp_fillnulls_14798> does not exist
Rename raster <usermask_mask.14798> to
[Raster MASK present]

On Thu, Dec 1, 2011 at 1:56 AM, Mark Seibel <mseibel@gmail.com> wrote:

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.

Maybe apply a raster mask over the ocean area?

Mark

Pankaj_Kr_Sharma · December 3, 2011, 7:39am

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.

--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/r-fillnulls-for-large-area-tp7047337p7057352.html
Sent from the Grass - Users mailing list archive at Nabble.com.

Nick_Jachowski · December 25, 2011, 4:12am

Dear All,

Thanks for all the help on this matter. Sorry for the late reply (fieldwork and vacation got in the way…)!

I ended up solving the problem by (1) filling the ocean nulls, (2) splitting the region into subregions, and (3) running it on a more powerful machine.

Best,
Nick

On Fri, Dec 2, 2011 at 11:39 PM, pkscwc <pkscwc@gmail.com> wrote:

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.

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:

If the system lacks RAM, then increase patience.

If the system lacks hard-disk space, don’t go for GIS.

Learn python.

–
View this message in context: http://osgeo-org.1803224.n2.nabble.com/r-fillnulls-for-large-area-tp7047337p7057352.html
Sent from the Grass - Users mailing list archive at Nabble.com.

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