[GRASS-user] G_calloc error using r.in.xyz

Hamish hamish_b at yahoo.com
Sat Feb 11 23:50:58 EST 2012

Patrick wrote:
> I've been able to get the attached visualizations of the attached data,

note these mailing lists go out to maybe 1000+ people, so there are strict
file size limits to attachments. Due to limitations in tcl/tk we still
haven't worked around NVIZ only exports as PPM or TIFF files. convert to
something 20+ times smaller on UNIXy computers with:
   convert image.tif new_image.png

(requires the ImageMagick package)

for larger images please consider posting a URL to a local server or
a service such as flickr, or for large data files something like

> however I'm still not quite sure how to use statistical tools to smooth
> it.

I don't think you need to smooth it.

Perhaps you set nulls to 0 (in the method=mean map), instead of 0s (in the
method=n map) to nulls?  i.e. it's not very noisy data, it's that instead
of skipping over null cells it's stretching the surface to the floor every
time it sees one.

as a comparison, this is what I get with your sample data:


# scan lidar xyz file for data extent:
GRASS> r.in.xyz -g in=pts_nocars_cut2.txt out=dummy fs=,
n=4133900.000000 s=3944500.000000 e=2714000.000000 w=2299200.000000 b=157340.000000 t=170330.000000

# data seems to be in intervals of 100 map units, so use that as our
# first guess as a cell resolution

# apply to region
GRASS> g.region -p -a res=100 n=4133900 s=3944500 e=2714000 w=2299200
north:      4133900
south:      3944500
west:       2299200
east:       2714000
nsres:      100
ewres:      100
rows:       1894
cols:       4148
cells:      7856312

# 1900x4150 is not too huge for an array size, so now let's try an import
# using method=n and see what the coverage is like at that resolution

GRASS> r.in.xyz in=pts_nocars_cut2.txt out=nocars.100.n method=n fs=,
r.in.xyz complete. 264876 points found in region.

# ignore cells with no hits (convert all 0s to NULLs)
GRASS> r.null nocars.100.n setnull=0

# check the stats:
GRASS> r.univar nocars.100.n
minimum: 1
maximum: 2
range: 1
mean: 1.06069
standard deviation: 0.238757

ok so hardly any cells have 2 hits in them, and no cells have more
than two cells in them. you've got a choice: keep what you've got and
interpolate between the dots to make a surface, or lower the cell
resolution and have less interpolation (basically educated guessing)
about what's in the gaps.

# let's try 500 map unit resolution and see what the "n" stats look like:
GRASS> g.region res=500 -a
GRASS> r.in.xyz in=pts_nocars_cut2.txt out=nocars.500.n method=n fs=,
GRASS> r.null nocars.500.n setnull=0
GRASS> r.univar nocars.500.n
n: 60060
minimum: 1
maximum: 10
range: 9
mean: 4.41019
standard deviation: 1.60499

ok, much better, an average of 4.4 hits per cell. the amount of hits
per cell will act as your smoothing filter on the raw data. (see r.in.xyz's various methods for more advanced versions than just taking the average)

if you need to satisfy statistical power you might need to have a higher
"n". for basic visualization let your eye decide, but as long as the
data isn't crap a couple or a few hits per cell should ensure that most
of the cells are not empty space between the hits.

happy that the grid array setup meets our needs, we can now proceed to
the actual data import:

GRASS> r.in.xyz in=pts_nocars_cut2.txt out=nocars.500.mean method=mean fs=,

In viewing the result we see the flight line is mostly covered, just the
occasional cell with no data in it. for visualization we want to fill the
tiny holes but not interpolate over places the plane didn't fly.

# first make a buffer with distance of 1 cell (r.grow would be fine too)

GRASS> r.buffer in=nocars.500.mean out=nocars.500.buffer dist=500

display it, make sure there are no empty holes where you want data
(red and yellow are ok areas). if there are lots of gaps you'd like
filled increase the buffer (or grow) distance a little. for the nocars
test case a buffer distance of 500 map units is fine, and we'll just
keep this map for later.

# now we fill the nulls:

GRASS> r.fillnulls in=nocars.500.mean out=nocars.500.filled

this can take a little while, be patient. if it takes more than a
lunchtime to complete consider alternate methods.. e.g. try the
r.surf.nnbathy module from grass addons or pester the grass developers
to backport the alternate v.surf.bspline method into grass 6.x.
A quick and dirty way is with r.surf.idw (or r.surf.idw2), but the
quality won't be as good as one of the spine fits.

once the hole filling is done you'll want to apply the valid-area MASK,

GRASS> r.mask in=nocars.500.buffer

# save a new map with the out-of-coverage areas cropped away:

GRASS> r.mapcalc "nocars.500.cropped = nocars.500.filled"

# now you can remove the MASK,

GRASS> r.mask -r

now you should have something pretty to look at,

GRASS> nviz nocars.500.cropped

go to the Visualize -> Raster Surfaces settings and change the Fine
resolution to 1 to show every cell. (depending on how huge your region
is, or if you want to apply any smoothing, you might make that show
every e.g. 3rd cell)

> Also, our ultimate viewing angle will be at low elevation and looking
> down the road, as if one were driving; a great deliverable would be an
> animation of someone moving along this roadway as if from a vehicles
> point of view.

that shouldn't be too hard, see  http://grass.osgeo.org/wiki/Movies

(or use the new & improved m.nviz.image command to make a series of static
images, then use an external encoder [e.g. mencoder or gifsicle] to make
the movie for a series of PNG input images)

> Is there a way to make it so that the points nearest
> you are well defined and as they get further away from your viewing
> angle they become statistically averaged/resolution decreases so that
> they look good and are there but are not eating up computer?

not really, or that is to say, probably OpenGL and your graphics card
take care of that a little bit for you. In practice NVIZ is really good
and it's not too much of an issue.

good luck,

More information about the grass-user mailing list