[GRASS-user] BIN raster import

[Glynn Clements]

Dylan Beaudette wrote:

> > note that GRASS considers the region bound coordinate to be at the
> > outer edge of the border cells.
> 
> Now that you mention this, is there an authoritative description of how GRASS 
> (and maybe GDAL) treat cells/pixels:

Inevitably, the authoritative description is the source code. Any
other documentation describes how modules are supposed to behave,
while the source code describes how they actually behave.

> 1. region calculations (outer edge of the border cells)

Well, the region isn't limited to raster data; it may also affect some
vector operations.

The region's bounds describe a rectangle in two-dimensional space. For
raster operations, this rectangle is subdivided into a grid of
rectangular cells, so the region's bounds are aligned with the edges
of the outermost cells.

> 2. cell locations (??? center, edge ???)

Cells are areas, not points, so they don't have locations. Their
corners have locations, as do their centres.

A cell with array indices (i,j) (easting, northing) corresponds to the
rectangle:

	{ (x,y) : west + i * ewres <= x < west + (i+1) * ewres,
		  north - (j+1) * nsres <= y < north - j * nsres }

whose centre is at:

	(west + (i+1/2) * ewres, north - (j+1/2) * nsres)

[Subject to wrapping of longitude values in lat/lon locations.]

> 3. raster to vector conversions (??? center, edge ???)

IIRC, r.to.vect uses the midpoints of the cell's edges (i.e. one
coordinate will be on a grid line, the other will be mid-way between
grid lines).

> 4. resampling (??? center, edge ???)

The built-in nearest-neighbour resampling of raster data calculates
the centre of each region cell, and takes the value of the raster cell
in which that point falls.

If the point falls exactly upon a grid line, the exact result will be
determined by the direction of any rounding error.

[One consequence of this is that downsampling by a factor which is an
even integer will always sample exactly on the boundary between cells,
meaning that the result is ill-defined.]

r.resample uses the built-in resampling, so it should produce
identical results.

r.resamp.interp method=nearest uses the same algorithm, but not the
same code, so it may not produce identical results in cases which are
decided by the rounding of floating-point numbers.

For method=bilinear and method=bicubic, the raster values are treated
as samples at each raster cell's centre, defining a piecewise-
continuous surface. The resulting raster values are obtained by
sampling the surface at each region cell's centre.

As the algorithm only interpolates, and doesn't extrapolate, a margin
of 0.5 (for bilinear) or 1.5 (for bicubic) cells is lost from the
extent of the original raster. Any samples taken within this margin
will be null.

AFAIK, r.resamp.rst behaves similarly, i.e. it computes a surface
assuming that the values are samples at each raster cell's centre, and
samples the surface at each region cell's centre.

For r.resamp.stats without -w, the value of each region cell is the
chosen aggregate of the values from all of the raster cells whose
centres fall within the bounds of the region cell.

With -w, the samples are weighted according to the proportion of the
raster cell which falls within the bounds of the region cell, so the
result is normally[1] unaffected by rounding error (a miniscule
difference in the position of the boundary results in the addition or
subtraction of a sample weighted by a miniscule factor).

[1] The min and max aggregates can't use weights, so -w has no effect
for those.

> I have often second-guessed myself on these very topics...

For the most part, the interpretation is the "obvious" one, given:

1. Cells are areas rather than points.

2. Operations which need a point (e.g. interpolation) use the cell's
centre.

>From a programming perspective, the functions:

	G_row_to_northing()
	G_col_to_easting()
	G_northing_to_row()
	G_easting_to_col()

all transform floating-point values.

Passing integer row or column indices to the first two functions will
return the coordinates of the cell's top-left corner; for the centre
coordinates, pass row+0.5 and/or col+0.5.

Similarly, the last two functions will typically return non-integral
values; use floor() to discard the fractional part and obtain the row
or column index of the cell within which the point lies.

-- 
Glynn Clements <glynn at gclements.plus.com>