[GRASS-dev] v.transform, wishlist #420

[Hamish]

Harri Kiiskinen wrote:
> That's what I thought when I looked at that code, but the thing is not
> so simple, as when using v.transform, there is not way to know,
> whether the 'to' or 'from' are coordinates are lat/long or not,
> because v.transform is used to transform between different coordinate
> systems, and the current location does not necessarily reflect either
> of these.

No, v.transform Is Not for transforming between different coordinate
systems. Use v.proj for that. It's for adjusting coordinates in 3D
space.

It's description implies otherwise: "Transforms a vector map layer from
one coordinate system into another coordinate system." but it's a bad
bad bad habit to get into and corrupts the whole grass location
paradigm. IMO the description should be changed. The r.region module
description is better: "Sets the boundary definitions for a raster map."
No more, no less.

v.transform and r.region are for low-level adjusting of a map's
coordinates. All maps within the same location are expected to use the
same coordinate space, without exception.

In practice v.transform may be used as a hack for transforming between
different coordinate systems, but really it would be better to import
your local data into a simple xy or custom projection, then pull/push
the data into your target location with v.proj and r.proj.

The problem is that a POINTS file with GCPs defines the relationship to
a projection, not a projection itself.

Currently for raster maps we have the i.rectify module which pushes
raster maps into another location; r.proj and v.proj which pull from
another location; the tcl gui georect tool which pulls maps from another
location (hopefully from a XY location); and v.transform which tries to
reproject in-place. This is a big pile of broken mess that needs to be
sorted out IMHO, and I'm pretty sure that the in-place solution is the
worst of the choices.

->We need to present the user with a single consistent process for
rectification/reprojection.


> An example:
> 
> I have a location with a local coordinate system ( x = [6000,9000], y
> = [2000,5000] ) for use with an archaeological excavation. To 'fix'
> this system to a global system, some points are measured also in
> lat/long coordinates. So I have a group of points having coordinates
> in the local, metric system and in the global, lat/long system. Just
> what v.transform needs to perform the transformation.

Be careful as it only works if your study area is very small; as lines
of longitude are not parallel.

> However, if I do the transformation in the original location, where
> the vector map is, the lat/long coordinates never go to G_lat_scan, as
> the 'projection' is not PROJECTION_LL. On the other hand, if I export
> the map and import it to the result location and do the transform
> there, the metric coordinates will be given to G_lat_scan(), as the
> 'projection' is PROJECTION_LL, even though for these points it is not
> true, and the resulting coordinates will suffer the normalization to
> the ranges mentioned above. I did test this already.

The easiest solution is to convert your lat/lon GCPs into your local UTM
zone coords (or other regional meter based system), then load your
custom system data into the UTM location, and then use v.transform.
Lat/lon is a pain for fine scale stuff anyway.

Perhaps ask on the PROJ4 mailing list for help on better defining a
custom projection if you wish. It's a fairly simple and empowering
exercise. Then you could use v.proj.

> So I still think, that there needs to be an additional test, either
> defined by the user or based on the probability, that ':' in the
> coordinate string means a lat/long that needs to be converted to
> decimal.

Applying radial (degrees) change in a cartesian coordinate system is not
correct! Euclidean space != polar space. Adding hacks to allow the
mixing of the two (only possible because the computer doesn't know
better) is a path to pain.

> > could the derived transform parms be used as the +towgs84 7-term
> > projection in the new location, then v.proj?
> 
> After a quick look at lib/vector/transform, I'd say that not
> self-evidently. Someone with good enough knowledge of linear algebra
> might be able to figure out, what the different transformation and
> rotation coefficients actually mean, and if they could be used in
> +towgs84. As the TODO file in that directory says, the "means of
> generating the coefficients seems rather odd".

The means aren't so important, as long as they work well enough. Once
you have the 'ends' (coefficients), you have what you need to calculate
the 7-term +towgs84 parm. (maybe?)

+towgs84=
3x translation parameters (Dx, Dy, Dz)
3x rotational parameters  (Rx, Ry, Rz)
1x scalar factor
    
v.transform: 
xshift,yshift,zshift   Shifting value for xyz coordinates
xscale,yscale,zscale   Scaling factor for xyz coordinates
zrot                   Rotation around z axis in degrees CCW

so Dx, Dy, Dz is the same, you "just" need to transform [Rx, Ry, Rz,
scale] to [scaleX, scaleY, scaleZ, curlZ].   (is it possible??)

After that the only frustrating bit is that you may have to adjust the
+/- signs to get the reference frame convention right, but that's about
it.

see http://proj.maptools.org/gen_parms.html#towgs84


> Then again, v.transform is not really about re-projecting the data,
> only about changing the coordinates into a different system

It shifts them within the same system. It is quite deceptive. The libgis
fn is not lacking, it's the module which is wrong in its presentation.

> (I think there is a difference?)

Reprojection *is* the changing of coordinates into a different system.
(and by system I mean precisely the +proj= param in the EPSG file, but
more generally as well, between two commonly used "systems" defined by
the entire EPSG code definition)


Note reprojection is not the same thing as rectification.

x1,y1 can be reprojected into x2,y2 if the projection systems for both
locations is known.  This is v.proj,r.proj = Projection->Projection.

Rectification (usually of remote imagery) consists of a perhaps uneven
or warped input image and ground control points. The GCPs are visually
applied by the user, and a mathematical fit between the two is derived
and applied, so the output then exists in a known projection. This is
i.rectify = XY location -> Projection. The mathematical relation is
found and then forgotten when the i.rectify module exits.
v.transform proably falls in this category (the mathematical relation is
found by the user and then not useful for much after v.transform is done).


> - usable just for the things its description says, putting dxf
> drawings in place, or, in my case, changing vectors from a very local,
> small scale coordinate system to something else - maybe another local
> system - just by defining some common points. I think it is nice just
> because it is simple, but I wouldn't use it for geographically large
> areas.

Right. Again a good solution is to transform your lat/lon points -> UTM
or similar, then use v.transform in the UTM location.

I think it is better to modify your method so it works within the
existing GRASS paradigm, and fix the way v.transform is presented, than
to muddy the GRASS paradigm even further by changing the way the libgis
fns and/or how modules work.


Sorry for the lengthy reply. If I were more of an expert in the ways and
had a real solution I'm sure this email would be much shorter. But you
raise an important problem that GRASS needs to address in a better and
more consistent way.


Hamish