Grid registration

From GRASS-Wiki
Jump to: navigation, search

Grid vs. cell-center registration

a subtle beast....


If you have gridded data -and do not need to interpolate- then r.in.xyz is the perfect tool, you just have to remember that GRASS uses cell-center registration not grid line registration.

  • if there are holes to interpolate you might as well just import as points directly to skip a r.to.vect step (assuming v.surf.rst)


So we want the height (depth?) points to be in the center of each cell in our region, but r.in.xyz's scan gives us the outer bounding box (ie outer grid lines).

Thus for gridded data you have to grow the scanned max extent by 1/2 a cell outwards if you want to use the max extent numbers.

A problem then is that 'g.region -a' and d.zoom will align the new region to the resolution & undoing your 1/2 cell. I'm not sure if the GUI zooms do that or not. After it gets damaged by d.zoom it renders as if the vector points are on the grid points, but really they aren't and halving the resolution or growing the region out by 1/2 a cell again gets them back to normal.

consider:

 n=1025  s=275  w=125  e=775  res=50

d.zoom (etc) will see the res=50 and want to grow it outwards to the next even multiple of res, so after d.zoom the region might be like:

 n=800  s=550  w=300  e=600   res=50

and the raster appears shifted by 1/2 a cell when you plot it.


In most cases this would be what you want to happen, but not when you have a grid-registration mismatch to deal with.

A q&d solution to all this is to halve the resolution again after import. BUT- do set the region back to the real value before starting NVIZ or else the results look really blocky. (or setting the fine resolution to "2" will do pretty much the same thing)


Another idea is to use r.region to shift the data by half a cell up and to the right. in this case the coordinate system is relative so we can set our "false eastings" and "y_0"s to anything we like.. 0,0 has little meaning in this coord system..

  • For something like sparse LIDAR data it doesn't really matter as there is no grid aliasing effect to worry about and r.in.xyz includes points falling exactly on the outer boundary of the region*
[*] except for the bottom boundary; I haven't gotten to that yet (see trac)


Example

# scan a text file containing gridded x,y,z data:
r.in.xyz -s -g in=tooth_scan.txt fs=, out=dummy
n=10.550000 s=3.850000 e=5.800000 w=1.050000 b=2.900000 t=4.875000


# apply what we know:
RES=0.05

# calc what half of that is:
RES2=`echo "$RES" | awk '{print $1 / 2}'`

# set the region based on that:
g.region n=10.550000 s=3.850000 e=5.800000 w=1.050000 res=$RES -p

north:      10.55
south:      3.85
west:       1.05
east:       5.8
nsres:      0.05
ewres:      0.05
rows:       134
cols:       95
cells:      12730


g.region n=n+$RES2 s=s-$RES2 w=w-$RES2 e=e+$RES2 -p

north:      10.575
south:      3.825
west:       1.025
east:       5.825
nsres:      0.05
ewres:      0.05
rows:       135
cols:       96
cells:      12960

# Now we have to be very careful, as 'g.region -a' and d.zoom
# will automatically align back to the grid. Save for backup:
g.region save=tooth

To check if the region alignment is set correctly, create a n map and enusre that you are in fact importing one point per cell:

r.in.xyz in=tooth_scan.txt fs=, out=tooth.n method=n --verbose
r.univar tooth.n

# if ''n'' is 1 thoughout, min,max,median, and mean maps should all be identical.
r.in.xyz in=tooth_scan.txt fs=, out=tooth method=mean

# the old looks-good-with-everything standby color table:
r.colors tooth color=bcyr
# my new favorite color table (in svn addons)
r.colors tooth col=haxby -n
# something more natural:
r.colors tooth color=sepia

# 'r.colors -e' can help to bring out the pic a bit

d.grid 1  # defaults are just right here


####

v.in.ascii in=tooth_scan.txt fs=, out=Tooth x=1 y=2 z=3

d.vect Tooth -z size=2 zcolor=bcyr

g.region res=$RES2 -p

d.zoom ...

####

# first thing after setting the z-exag to 1.0 go into the lighting
#  controls and set lights to follow viewpoint or else the shadows getcha
g.region tooth
nviz tooth