Subscribe to this thread
Home - General / All posts - Computing <geom> in epsg::3857
joebocop
348 post(s)
#07-Nov-18 23:15

I have a table, [Point of Interest], which has epsg::4326 coordinate columns [latitude_wgs84] and [longitude_wgs84].

I would like to add two computed <geom> type fields to the table; one for epsg::4326 geometry, and the other for epsg::3857 geometry.

I have found that the following SQL succeeds, though I suspect I am doing something wrong.

--$manifold

ALTER TABLE [Point Of Interest]

(

  ADD [geom_ll_mfd] GEOM AS [[ 

    GeomMakePoint(VectorMakeX2([longitude_wgs84],[latitude_wgs84])) 

  ]],

  ,ADD INDEX [geom_ll_mfd_x] RTREE ([geom_ll_mfd])

  ,ADD PROPERTY 'FieldCoordSystem.geom_ll_mfd' CoordSystem(4208)

  ,ADD [geom_wm_mfd] GEOM AS [[ 

    CoordConvert(

      CALL CoordConverterMake(CoordSystem(6505), CoordSystem(4208)), 

      GeomMakePoint(VectorMakeX2([latitude_wgs84][longitude_wgs84]))

    ) 

  ]]

  ,ADD INDEX [geom_wm_mfd_x] RTREE ([geom_wm_mfd])

  ,ADD PROPERTY 'FieldCoordSystem.geom_wm_mfd' CoordSystem(6505)

);

Why do I have to specify the coordinates in opposite order for the two GeomMakePoint() functions?

The CoordConverterMake() function contains the definitions for the source and target coordinate systems, so shouldn't it assume the "source" geom's coordsys is epsg::4326 (Manifold 4208), having Y/X coordinate order?

Again, I suspect I've done something wrong and have accidentally suceeded here, so please don't hesistate to correct me. Thanks.

tjhb

8,335 post(s)
#08-Nov-18 01:18

CoordConverterMake() takes arguments in the order <target>, <source>. (I always have to remind myself of this too.)*

Regardless, though, are you sure you want to do this with a computed field? It's inefficient to create a new CoordConverter object for every row (then throw it away). A static query could reuse the object, making it much more efficient. But of course, that would not reproject on-the-fly.

*[Added.] That is one thing. The other is this that EPSG::4326 does indeed use YX coordinate order.

joebocop
348 post(s)
#08-Nov-18 04:16

A static query is certainly how I have this same behaviour wired up in release 8, but am I wrong in believing that release 9 favours computed fields over queries returning new geoms?

It's inconvenient to manually generate temporary spatial index on drawings from queries, so I'd been favouring computed geom columns, or even static geom columns regenerated by an UPDATE query on demand. Maybe I need to look at Drawings from queries again.

edit: There's also no ability to edit data in fields returned by a query in cases such as these, where new fields are created by the query.

adamw


8,204 post(s)
#08-Nov-18 15:00

The CoordConverterMake() function contains the definitions for the source and target coordinate systems, so shouldn't it assume the "source" geom's coordsys is epsg::4326 (Manifold 4208), having Y/X coordinate order?

This is exactly what happens - you say that source coordinate system is EPSG:4326, which uses YX coordinate order, so when you specify source geom, you have to pass Y (lat) first and X (lon) second.

Good code, by the way. Two notes:

1. If what you want the query to do is convert between EPSG coordinate systems, then instead of doing CoordSystem(xxx) and remembering that 4208 is EPSG:4326 and 6506 is EPSG:3857, you can just write 'EPSG:4326' or 'EPSG:3857', eg:

--SQL9

...

CoordConvert(

  CALL CoordConverterMake('EPSG:3857''EPSG:4326'),

  GeomMakePoint(VectorMakeX2(...))

)

...

2. You can avoid creating a coordinate converter object every time the system has to recompute the value of 'geom_wm_mfd' by using a shared object:

--SQL9

 

CREATE TABLE [poi] (

  [mfd_id] INT64,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

  [longitude_wgs84] FLOAT64,

  [latitude_wgs84] FLOAT64

);

 

INSERT INTO [poi] ([longitude_wgs84][latitude_wgs84])

VALUES (30, 10), (30, 20), (30, 30), (30, 40), (30, 50);

 

ALTER TABLE [poi] (

  ADD [geom_ll_mfd] FLOAT64X2 AS [[ 

    VectorMakeX2([longitude_wgs84][latitude_wgs84])

  ]],

  ADD [geom_wm_mfd] FLOAT64X2 WITH [[

    VALUE @conv TABLE = CALL CoordConverterMake('EPSG:3857''EPSG:4326');

  ]] AS [[

    CoordConvertPoint(@conv,

      VectorMakeX2([latitude_wgs84][longitude_wgs84]))

  ]]

);

I used FLOAT64X2 instead of GEOM to let you see the values directly in the table window.

adamw


8,204 post(s)
#08-Nov-18 15:05

Forgot to add:

If the coordinate system of 'geom_ll_mfd' is supposed to be EPSG:4326, you actually have to reverse lat/lon values there as well. I did not do this in the example above, so in my example - and in your query - the coordinate system of 'geom_ll_mfd' is not EPSG:4326, but rather plain 'Latitude / Longitude' which uses XY coordinate order.

Manifold User Community Use Agreement Copyright (C) 2007-2017 Manifold Software Limited. All rights reserved.