g.setproj, r.proj, rant about s.surf.***

Simon Cox s.cox at dem.csiro.au
Sat Aug 27 09:07:23 EDT 1994


>
>I have 4.1 but I don't (generally) do vectors.  Is there an easy way
>in GRASS to convert raster data between different projections?  Really
>I'm just curious as I would like to be able to overlay images from
>different sources with different projections.  I'm thinking of adding
>this functionality into Khoros v2, but would be interested in knowing
>if it already exists in GRASS or elsewhere.
>
>Tim -- "If only the world were flat!"

Ah - this is the (by now getting old) r.proj problem.  To my knowledge this
has not yet been written in a robust enough form to be included in any of
the release.  Last time I needed to do this I think I did something like
...

(i)     starting in GRASS in the location that the original raster is in,
write it out as a site-list using
r.stats -1gz ...

(ii)    exit from GRASS.  Use "proj" to convert the site list into the
projection that you need (you may have to have an intermediate step of
lat-long).  You may need also to use "awk" to reformat the columns at some
stage.

(iii)   fire up GRASS in the new location/projection etc.

(iv)    import the site-list using s.in.ascii.

(v)     interpolate the site-list onto a raster using s.surf.idw or
s.surf.tps.

For this last step, note that (a) due to its lack of intelligent indexing
of the site list, s.surf.idw is excruciatingly slow for big data sets (more
than a few hundred);  (b) however, since you _will_ want to produce a final
raster at a similar resolution to the original one, you don't need to use
many points to generate the new values - ie you can set npmin in s.surf.tps
to a very small number (eg 10) to get a good result.

However, I still think that using s.surf.tps is silly for a task like this,
where a much simpler (eg bilinear) interpolation would be more than
adequate - the new cell positions are surrounded by "sites" from the
original cells.  The only reason I used tps was because it actually ran
_faster_ than idw, due to the fact that an indexing scheme is used properly
for tps!  I understand why the effort has gone into these clever
interpolators, and am thankful that tps is available.  But I think the fact
that there is not a fast bilinear or nearest neighbour s.surf.*** program
is a serious deficiency in GRASS's armory.  At the very least, couldn't
s.surf.idw incorpporate the quadtree from tps, or maybe a more general
s.surf program could be written with -method=tps|idw|nn|bilin as an option,
so that large chunks of the code are reused?

(I'm sorry I don't have time to just do it myself.)

When this has all been done, it should be possible to encapsulate the above
in a script ...

End of rant.            Simon Cox

____
__________________________________________________________________
Dr Simon Cox                             __  L
CSIRO Exploration & Mining            ,~'  L_|\       Australian
39 Fairway, PO Box 437,            ;-'         \      Geodynamics
Nedlands, WA  6009,  Australia     (            \     Cooperative
      Phone +61 9 389 8421         +    ___     /     Research
      Fax   +61 9 389 1906          L~~'   "\__/      Centre
s.cox at dem.csiro.au                            7
__________________________________________________________________





More information about the grass-user mailing list