Help with understanding the localframe method

[Tom V]

Hello

I do not understand the Metashape.CoordinateSystem.localframe() method. This may well be because I have a poor understanding of the underlying coordinate systems being used. If anybody could shine any light on this I would be really grateful!

The manual tells us that localframe() "Returns 4x4 transformation matrix to LSE coordinates at the given point". The input parameter is a vector that is "Coordinates of the origin in the geocentric coordinates" and it then returns a transformation matrix that transforms "from geocentric coordinates to local coordinates".

So I think I am probably misunderstanding the above definition. By my (incorrect) understanding, in the following code the region center at the end should be at the origin in the local coordinates i.e. (0,0,0).

Instead it is at ~ (0, 0, 2187)

Code: [Select]

import Metashape
doc = Metashape.app.document
chunk = doc.chunk
M = chunk.transform.matrix
world_center = M.mulp(chunk.region.center) # Convert region centre to a geocentric coordinate
T_origin = chunk.crs.localframe(world_center) # Produce a matrix (T_origin) that tranforms geocentric coordinates to local coordinate using the region center as the origin.

print(f"\nRegion centre in local coordinates before before transformation: {chunk.region.center}\n")
print(f"Region centre as world coordinate: {world_center}")

# Apply matrix to the region centre in geocentric coordinates
# By the definition of T_origin this should therefore be at the origin of the local coordinate system - but it is not!
chunk.region.center = T_origin.mulp(world_center) 

print(f"\nRegion centre in local coordinates after transformation: {chunk.region.center}")
Output:

>>> Region centre in local coordinates before before transformation: Vector([0.28900891187347366, 1.163145900441732, -10.650618341058543])
 
>>> Region centre as world coordinate: Vector([4244594.686245037, 959551.9924872194, 4649362.857761447])
 
>>> Region centre in local coordinates after transformation: Vector([2.546585164964199e-11, 0.0, 2187.11612998927])

Confused... 

Many thanks in advance!

Tom

[Paulo]

Hello Tom,

I would look at what is your chunk.crs? Noirmally it  should be WGS84 or some projected CS from WGS84 as you can see on following output, the LSE will have its X and Y tangent to geographic CS or PCS with X pointing East and Y north and Z will be equal to elipsoid height at this point.

Code: [Select]

crs
Out[10]: 2023-12-08 17:00:37 <CoordinateSystem 'WGS 84 (EPSG::4326)'>

O_geoc = ps.Vector([4244594.686245037, 959551.9924872194, 4649362.857761447])

O_geog = crs.project(O_geoc)

O_geog
Out[9]: 2023-12-08 16:57:42 Vector([12.738420342965075, 47.08595860158432, 1482.88178207708])

crs.localframe(O_geoc).mulp(O_geoc)
Out[13]: 2023-12-08 17:02:57 Vector([-2.3283064365386963e-10, 0.0, 1482.8817820767872])

[Tom V]

Thanks very much Paulo. Very much appreciated. I can follow your code but I think I am misunderstanding some basic underlying principles. I have spent ages reading around these subjects but I am clearly misunderstanding with respect to Metashape principles. I am sure this is explained somewhere but I have not found it. I am going to have to ask some really dumb questions...

So, my understanding is that in a typical georeferenced project using a projected coordinate system (EPSG::31255 in my case) Metashape deals with three systems of coordinates.

1. The internal coordinate system. This has arbitrary scale and orientation with respect to the real world and is defined by three orthogonal axes. Is this what is meant by "LSE coordinates"?

2. A geocentric coordinate system. This is a coordinate system based around a spheroid that is defined by the geographic coordinate system applicable to the chunk. In my case it is the Bessel 1841 (EPSG:7004) spheroid that is used with the projected coordinate system EPSG::31255. Conversion from the internal coordinates of (1) to this geocentric system is achived using the chunk.transform.matrix.

This geocentric system confuses me. What units are the values returned by this system? I would expect angular values if this is a sphere but they are values in the thousands, which doesn't seem to fit. I cannot relate these to anything real world.

3. A projected coordinate system. Achieved by applying the crs.project() method to the geocentric coordinate system of (2). With this I am on safe ground  as it is either a projected coordinate in metres, relating to the crs of ths chunk, or an angular position in degrees if a geographic coordinate systemn is used. Either way it is familiar.

An example of this confusion. When I apply the matrix from the localframe method to a geocentric coordinate (2) I am expecting a new vector that is in the internal coordinate system (1). Yet, the scaling of the locaframe matrix is 1 (i.e. no scaling) whereas there is a scaling > 10 involved in the chunk.transform.matrix, which also converts between these two coordinate systems. If it scales converting one way then it should inversely scale going the other way. So, I must have something wrong.

These coordinate system relationships have confused me for ages and any thoughts would be hugely appreciated!

Thanks in advance!

Tom

[Paulo]

Hi,

no he internal coordinate system is not the LSE CS. For a given point O, the LSE CS is tangent  to the ellipsoid at position O with Xlse pointing east, Ylse point north and Zlse up. However it is is tangent at ellipsoid with h or Z 0..

In example the blue axes are the LSE CS for point region center but at 0 height while project is at h 268...the red axes are the internal coordinate system.

So if you have:

T = chunk.transform.matrix
L = crs.localframe(Ogeoc) where Ogeoc are the geocentric coordinates foe a given point O

then T goes from Internal CS to Geocentric CS and L goes from  Geocentric CS to LSE CS for given point O.

As for geocentric CS it is a cartesian orthogonal centered at earth center and values in meters with Xgeoc pointing to Greenwich meridian and Zgeoc to North Pole... That is why values are in the millions as on earth surface we are sone 6 million meters from center see 2nd attachment...

So if you transform geographic coordinates (long 0, lat 0, h 0 )) to geocentric you will get radius of elipsoid at equator for Xgeoc. And from lng 0, lat 90 and h 0 you will get radius at North pole (semi minor radius) for Zgeoc.

Code: [Select]

In [10]: crs.transform(ps.Vector((0,90,0)),wgs84,wgs84.geoccs)
Out[10]: 2023-12-09 17:06:13 Vector([9.52360149646502e-09, 0.0, 6356752.314245179])

In [11]: crs.transform(ps.Vector((0,0,0)),wgs84,wgs84.geoccs)
Out[11]: 2023-12-09 17:06:44 Vector([6378137.0, 0.0, 0.0])

Hope this helps...

[Tom V]

Paul, that was extremely helpful, thank you! Really clear.

So many things that were previously incoherent to me, both in this thread and before, now make sense!

Next is to understand the propagation of variance in relation to the covariance matrix (I have seen your post but need to read up on the maths)... but that's for another day!

Thanks so much 

Tom