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

GRASS GIS GRASS GIS developer
1

Hello all,

on advise from Maciek, I'll re-post here the comments I added to the
wishlist item mentioned in the Subject. The item is about v.transform
and how it does not accept lat/long in ascii format in the 'pointslist'
file. Hamish responded, that it probably would be easy to fix by just
modifying the module to use G_scan_easting() and G_scan_northing().
However, and here I quote myself:

<quote>
I looked into the matter somewhat deeper, and it turns out not
to be as simple as just putting in G_scan_*()-functions. These
functions delegate the work to ll_scan(), through G_lat_scan()
and G_long_scan(), and ll_scan() normalizes its resulting arguments
to the ranges [-180,180] and [-90,90]. The result is, that the
coordinates can be given to these functions only if it is supposed
to be lat/long-coordinate.

As I have some interest for this module, this might be a chance
for a first code contribution to the project. I have two possible
options for the solution:

1. Check the coordinate string first to see if the character
':' can be found in it, and if it can be, assume the coordinate
is lat/long and do the necessary conversion

2. Add two more flags to the module, to indicate that the from-
and/or to-coordinates should be considered as lat/long.

The first solution is more elegant, but the second one makes
sure, that the user knows the form of the coordinates.

Two questions:

Which of these would be the preferred solution?

Is is ok if I meddle with this module and send in diffs/patches
at some point?

In the long run, I'd like to further modify this module so, that the
resulting vector could be placed in a different location/mapset, so it
could be used to transform vector files from a location based on a
local coordinate system to another, 'global' system.
</quote>

Harri K.

2

Harri Kiiskinen wrote:

I looked into the matter somewhat deeper, and it turns out not
to be as simple as just putting in G_scan_*()-functions. These
functions delegate the work to ll_scan(), through G_lat_scan()
and G_long_scan(), and ll_scan() normalizes its resulting arguments
to the ranges [-180,180] and [-90,90]. The result is, that the
coordinates can be given to these functions only if it is supposed
to be lat/long-coordinate.

lib/gis/wind_scan.c

int G_scan_northing ( const char *buf, double *northing, int projection)
{
    if (projection == PROJECTION_LL)
    {
        if(G_lat_scan (buf, northing))
            return 1;
        if (!scan_double (buf, northing))
            return 0;
        return (*northing <= 90.0 && *northing >= -90.0);
    }

    return scan_double (buf, northing);
}

it only calls G_lat_scan() if the location's projection is lat/lon.
No need for any extra checks.

In the long run, I'd like to further modify this module so, that the
resulting vector could be placed in a different location/mapset, so it
could be used to transform vector files from a location based on a
local coordinate system to another, 'global' system.

could the derived transform parms be used as the +towgs84 7-term
projection in the new location, then v.proj?

Hamish

3

On Tue, 2007-06-12 at 16:37 +1200, Hamish wrote:

Harri Kiiskinen wrote:
> I looked into the matter somewhat deeper, and it turns out not
> to be as simple as just putting in G_scan_*()-functions. These
> functions delegate the work to ll_scan(), through G_lat_scan()
> and G_long_scan(), and ll_scan() normalizes its resulting arguments
> to the ranges [-180,180] and [-90,90]. The result is, that the
> coordinates can be given to these functions only if it is supposed
> to be lat/long-coordinate.

lib/gis/wind_scan.c

int G_scan_northing ( const char *buf, double *northing, int projection)
{
    if (projection == PROJECTION_LL)
    {
        if(G_lat_scan (buf, northing))
            return 1;
        if (!scan_double (buf, northing))
            return 0;
        return (*northing <= 90.0 && *northing >= -90.0);
    }

    return scan_double (buf, northing);
}

it only calls G_lat_scan() if the location's projection is lat/lon.
No need for any extra checks.

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.
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.
  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.
  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.

> In the long run, I'd like to further modify this module so, that the
> resulting vector could be placed in a different location/mapset, so it
> could be used to transform vector files from a location based on a
> local coordinate system to another, 'global' system.

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".
  Then again, v.transform is not really about re-projecting the data,
only about changing the coordinates into a different system (I think
there is a difference?) - 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.

Harri K.

4

Hi Harri,

It sounds like you are interested in a projection-based
transformation, rather than the affine transformation which
v.transform can apply.

Here are the steps that I would take:

1. collect control points, using a projected coordinate system

custom coordinates
( x,y )

UTM for example
( x', y' )

2. apply affine transform in UTM coordinates

3. inverse project to geographic coordinates, if that is what you really want.

You will need:

1. a new location in UTM
2. import your custom coordinates into that location, ignoring the
fact that the coordinate systems do not match
3. use v.transform

I think that v.transform should issue some kind of error if invoked
within a LL location...

cheers,

Dylan

On 6/12/07, Harri Kiiskinen <harkiisk@utu.fi> wrote:

On Tue, 2007-06-12 at 16:37 +1200, Hamish wrote:
> Harri Kiiskinen wrote:
> > I looked into the matter somewhat deeper, and it turns out not
> > to be as simple as just putting in G_scan_*()-functions. These
> > functions delegate the work to ll_scan(), through G_lat_scan()
> > and G_long_scan(), and ll_scan() normalizes its resulting arguments
> > to the ranges [-180,180] and [-90,90]. The result is, that the
> > coordinates can be given to these functions only if it is supposed
> > to be lat/long-coordinate.
>
> lib/gis/wind_scan.c
>
> int G_scan_northing ( const char *buf, double *northing, int projection)
> {
> if (projection == PROJECTION_LL)
> {
> if(G_lat_scan (buf, northing))
> return 1;
> if (!scan_double (buf, northing))
> return 0;
> return (*northing <= 90.0 && *northing >= -90.0);
> }
>
> return scan_double (buf, northing);
> }
>
> it only calls G_lat_scan() if the location's projection is lat/lon.
> No need for any extra checks.

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.
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.
        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.
        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.

> > In the long run, I'd like to further modify this module so, that the
> > resulting vector could be placed in a different location/mapset, so it
> > could be used to transform vector files from a location based on a
> > local coordinate system to another, 'global' system.
>
> 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".
        Then again, v.transform is not really about re-projecting the data,
only about changing the coordinates into a different system (I think
there is a difference?) - 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.

Harri K.

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

5

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

6

Dylan Beaudette wrote:

You will need:

1. a new location in UTM
2. import your custom coordinates into that location, ignoring the
fact that the coordinate systems do not match
3. use v.transform

of course Dylan got to the same solution first.. :slight_smile:

I think that v.transform should issue some kind of error if invoked
within a LL location...

but, what if due to a transcription error or whatever you want to shift
everything a degree to the west?

(I do have a chart here where all lons are off by 1". It's a hand drawn
1853 nautical chart from the HMS Acheron which was only resurveyed and
replaced in the last ~5 years. I met a guy from the British Admiralty
service on a dive trip to the great barrier reef soon after- he said
they had a big party as it was the last of their old-old maps still in
service to be replaced. I guess that makes us the ends of the known
world down here :slight_smile:

Hamish

7

Hamish 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."

Whoever wrote that was presumably thinking about the use of the term
"coordinate system" in mathematics/geometry rather than in cartography.

As the user will probably view it from a cartographic perspective,
it's rather misleading.

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

8

Thanks for Glynn and Hamish for clarifications; I fully accept your
reasoning. I think I got slightly carried out by what could be done and
thus ignoring whether it really should be done.

Harri K.