Subscribe to this thread
Home - General / All posts - M9 tile addition and subtraction in surfaces
danb

2,064 post(s)
#28-Feb-19 03:32

I have two image components representing (in this case) a geoid model and an ellipsoidal height lidar DSM which have the following properties:

1. Name: [new_zealand_quasigeoid_2016_raster]

CoordSys: { "Accuracy": 1, "Axes": "YX", "Base": "NZGD2000 (EPSG:4167)", "Eccentricity": 0.08181919104281579, "LocalOffsetX": 159.99166665, "LocalOffsetY": -63.39167345, "LocalScaleX": 0.0166667, "LocalScaleY": 0.0166667, "MajorAxis": 6378137, "Name": "NZGD2000 (EPSG:4167)", "Rect": [ 160.6, -55.95, -171.2, -25.88 ], "System": "Latitude \/ Longitude", "Transform": "Molodensky-Badekas", "Unit": "Degree", "UnitLatLon": true, "UnitScale": 1, "UnitShort": "deg" }

TileSize: [ 256, 256 ]

Rect: [ 0, 203, 1801, 2304 ]

2. Name: [T3305840 Interpolate, Kriging Image]

CoordSys: { "Name": "Universal Transverse Mercator Zone 60 (S)", "System": "Transverse Mercator", "CenterLat": 0, "CenterLon": 177, "FalseEasting": 500000, "FalseNorthing": 10000000, "ScaleX": 0.9996, "ScaleY": 0.9996, "Axes": "XYH", "Base": "World Geodetic 1984 (WGS84)", "MajorAxis": 6378137, "Eccentricity": 0.08181919084262149, "Unit": "Meter", "UnitScale": 1, "UnitShort": "m" }

TileSize: [ 128, 128 ]

Rect: [ 331023, 5840000, 332000, 5842000 ]

So basically two images in different projections with different offsets, tile sizes and pixel dimensions.

I want to be able to add or subtract the value of the [new_zealand_quasigeoid_2016_raster] pixels from the [T3305840 Interpolate, Kriging Image] pixels (I guess) using something along the lines of:

[T3305840 Interpolate, Kriging Image].[Tile] - [new_zealand_quasigeoid_2016_raster].[Tile]

My feeling is that this process would be very fast within M9, but that I would first need to reproject and retile one or other of the images such that the tiles of both images were of the same size as each other and are projected into the same space so that I can join the images on tile X, Y.

If this indeed the way to approach it, I unfortunately have no idea of how to go about it so any guidance would be very much apperciated. I also have a sample project but this weights in at ~30mb but I can put this on an ftp if required.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

adamw


10,447 post(s)
#28-Feb-19 07:56

You can reproject images on the fly using a query.

Here is an example (I am not attaching test data, because even the MXB for what I used is over the 4 MB limit, if you have any difficulties in applying the query, tell me, and I will upload test data somewhere):

--SQL9

 

-- takes [busherh_reactor]

-- subtracts [busherh_reactor 2] (different tile size / projection)

-- writes [busherh_reactor 3] (same tile size / projection as first image)

--

-- all images are single-channel

 

VALUE @conv TABLE = CALL CoordConverterMake(

  ComponentCoordSystem([busherh_reactor]),  -- target, 128x128

  ComponentCoordSystem([busherh_reactor 2]-- source, anything

);

 

INSERT INTO [busherh_reactor 3] (x, y, tile)

SELECT x, y, [tile] -

  CoordConvertTile(@conv, [busherh_reactor 2],

    VectorMakeX4(x * 128, y * 128, (x+1) * 128, (y+1) * 128))

FROM [busherh_reactor];

 

TABLE CALL TileUpdatePyramids([busherh_reactor 3]);

The reprojection part is CoordConvertTile. You pass it a rect in the target coordinate system for which you want to compose a tile, a source image, and a coordinate converter, and the function composes a new tile with projected values. Tile size of the source image is irrelevant.

Now, I checked and we have an outstanding issue with CoordConvertTile in that on some combinations of coordinate systems it produces horizontal artifacts on tile edges - we'll re-check what's up with that and fix it.

danb

2,064 post(s)
#01-Mar-19 02:58

Many thanks for this Adam. Much simpler than I had anticipated and incredibly fast. I have some large grids which I will need to manipulate using variants of this method and this provides an extremely useful template.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

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