[Anuga-user] pixel vs gridline registration issues

[Ole N. <ole...@gm...>]

Hello Joaquim and other list subscribers

Thank you very much for the posting and for enlightening us on the topic of
image registration. To be totally honest, we have never thought about this
in our development of ANUGA.
Talking to the team, I am pretty convinced that we had the gridline paradig=
m
in mind when we wrote the code to export from triangular meshes to ASCII
grids. The format was inspired by ARC which is what we needed to export to
(and from).
Hence, we take the lower left corner (xllcorner, yllcorner) to be the
west-most and south-most point of the grid with all other points being
offset as multiples of the grid resolution:

x_i =3D xllcorner + i*offset
y_j =3D yllcorner + j*offset

with i, j =3D 0, 1, ..., N

This means that we would interpret the first points (i=3D0, j=3D0) as
  x_0 =3D xllcorner
  y_0 =3D yllcorner

We believe this is the way grids are interpreted by ARC but we could be
wrong. In any case it sounds like we are using grid line registration
according to the definitions you pointed to.

It is, of course, entirely possible that we have made mistakes in the
conversion as you suggest and we welcome  your  efforts to correct  these.
In fact, this was one of the motivations for making ANUGA open source: That
skilled people like yourself would pickup errors and inconsistencies so tha=
t
the model can improve over time.

If you were to rewrite the export routine, we would be very happy to includ=
e
it in ANUGA (wither appropriate acknowledgement, of course), but we would
ask that any new or improved routine comes with unit tests (as you can see
in test_data_manager.py).
Alternatively, the problem might be resolved by changing the easting_max
line you mention to follow the grid_line paradigm.

In any case, we are grateful for this kind of feedback  and eager to improv=
e
the ANUGA code base together with the growing community.

Best regards
Ole Nielsen
Geoscience Australia


----Original message
Hello,

greetings from the Atlantic.

First, a word of warning for anyone who might be using the GMT program
grd2xyz -E to convert into
Arc/Info .asc files. That program is bugged in the sense that it
produces wrong lower left corner coordinates.
It does create xllcenter & yllcenter keewords (which I'm afraid anuga
interprets as xllcorner & yllcorner)
but its coordinates are shifted 1/2 cellsize in the west and north
direction.

This issue leads me to the main point of this message the - pixel vs
gridline registration - confusion.
For a nice (with graphical) explanation of this matter see this GMT doc pag=
e
http://gmt.soest.hawaii.edu/gmt/doc/html/GMT_Docs/node158.html

By definition referencing with xllcenter & yllcenter means that the
pixel registration is adopted.
The convert_dem_from_ascii2netcdf() copies the .asc header into the .dem
netCDF file and so
we can say that the .dem grid is pixel registered.
Now the dem2pts() function. At line n=BA 1335 we have

if easting_max is None: easting_max =3D xllcorner + ncols*cellsize

this is correct in the pixel registration universe but I wonder if that
is really the intention
for at line 1372, we have now

for j in range(ncols):
x =3D j*cellsize + xllcorner

but these x are computed using the gridline registration paradigma.
Otherwise, it would have to be
x =3D (j + 0.5)*cellsize + xllcorner


Please correct me if I am wrong

Best regards

Joaquim Luis