Subscribe to this thread
Home - General / All posts - M9 Transfer Coordinate of Highest/Lowest Contained Pixel from Raster to Drawing
atrushwo59 post(s)
#31-Jul-19 18:13

Hi All,

I have a drawing containing areas [DRAWING] and a surface [GROUND]. I would like to transfer the coordinate of the lowest/highest elevation pixel within each area to the drawing. I had an (inelegant) way to do this in M8, but cannot figure it out in M9. My M8 SQL is below. Any ideas on how I might accomplish this in M9?

UPDATE

(

SELECT

 [1].[LP] AS [A],

 [POINT] AS [B]

FROM

(

SELECT

 [SOURCE],

 FIRST([POINT]AS [POINT]

FROM

 (

 SELECT 

 [1].[ID] AS [SOURCE],

 PROJECT(ASSIGNCOORDSYS(NEWPOINT([S].[CENTER X (I)],[S].[CENTER Y (I)]),COORDSYS("GROUND" AS COMPONENT)),COORDSYS("DRAWING" AS COMPONENT)) AS [POINT]

 FROM [DRAWING] AS [1]

 INNER JOIN [GROUND] AS [S] ON Touches([1].[ID],ASSIGNCOORDSYS(BoundingBox(NewLine(NewPoint([S].[X (I)][S].[Y (I)]),NewPoint([S].[X (I)] + 1, [S].[Y (I)] + 1))),COORDSYS("GROUND" AS COMPONENT)))

 ORDER BY [S].[HEIGHT (I)] ASC

 )

 GROUP BY [SOURCE]

)

INNER JOIN [DRAWING] AS [1] ON [SOURCE] = [1].[ID]

)

SET [A] = [B]

Dimitri


5,620 post(s)
#31-Jul-19 18:23

There's an illustrated, step by step topic in the user manual to show you how to do that. See the Example: Transfer DEM Terrain Heights to Areas in a Drawing topic

atrushwo59 post(s)
#31-Jul-19 18:28

Forgive me, but this appears to transfer the max/min elevations from the surface to the drawing. I want to transfer the coordinates of the max/min elevations from the surface to the drawing. I don't believe the example provided indicates how the latter is accomplished.

Dimitri


5,620 post(s)
#01-Aug-19 07:32

Read the topic carefully all the way to the end... it discusses transfer of min/max values. Getting the coordinate of that location is more difficult, but you can use a very similar approach to do that. Cover your areas with a set of smaller areas of desired resolution, do the same transfer to those, and then pick out the one within each bigger area with the highest min/max value.

tjhb

8,958 post(s)
#01-Aug-19 11:05

I think I have a better approach, much faster.

I will post the code tomorrow.

Of course, test data would be very helpful before then.

Dimitri


5,620 post(s)
#01-Aug-19 14:57

I'm looking forward to that. After thinking about it some more, what I suggested is a crude approach.

atrushwo59 post(s)
#01-Aug-19 15:13

Hey Tim,

I've attached a sample map. I think their could be a way to do it with a TileToValues command? Reconstructing the raster with a series of points seems rather crude though. I'm not sure there is a better alternative though?

Attachments:
map_TransferCoord.map

atrushwo59 post(s)
#01-Aug-19 20:48

I think I figured it out finally. Let me know if anyone has any improvements.

CREATE ROOT X;

USE ROOT X;

CREATE DATASOURCE D AS ROOT;

FUNCTION A_TBL() TABLE AS D::[WATERSHEDS] END;

FUNCTION A_DWG() TABLE AS D::[WATERSHEDS DRAWING] END;

FUNCTION S_TBL() TABLE AS D::[FILLED_A11 TILES] END;

FUNCTION S_DWG() TABLE AS D::[FILLED_A11] END;

VALUE @conv TABLE = CALL CoordConverterMake(

  ComponentCoordSystem(CALL A_DWG()),

  ComponentCoordSystem(CALL S_DWG()));

SELECT CoordConvert(@conv, GeomMakePoint(VectorMakeX2([x] * 128 + [x 2][y] * 128 + [y 2]))) AS [GEOM][VALUE] AS [HEIGHT]

INTO [TEMP]

FROM

(

SELECT [X],[Y],

SPLIT CALL TILETOVALUES([TILE])

FROM CALL S_TBL()

)

WHERE [VALUE] IS NOT NULL;

ALTER TABLE [TEMP] (ADD INDEX [GEOM_XX] RTREE ([GEOM]));

ALTER TABLE [TEMP] (ADD PROPERTY 'FieldCoordSystem.Geom' ComponentCoordSystem(CALL A_DWG()));

CREATE DRAWING [TEMP Drawing] (

  PROPERTY 'Table' '[TEMP]',

  PROPERTY 'FieldGeom' 'Geom'

);

SELECT

 [S_MFD_ID],

 [O_HEIGHT],

 [O_GEOM]

INTO [TEMP2]

FROM 

 CALL GeomOverlayContainedPar(

 CALL A_DWG() ([mfd_id],[geom]), --S

 [TEMP DRAWING] ([geom],[HEIGHT]), --O

 0,ThreadConfig(SystemCpuCount()))

ALTER TABLE [TEMP2] (ADD INDEX [S_XX] BTREEDUP ([S_MFD_ID]));

SELECT 

 [1].[S_MFD_ID],

 ((COLLECT [1].[O_GEOM] ORDER BY [1].[O_HEIGHT] ASC FETCH 1))

INTO [TEMP3]

FROM [TEMP2] AS [1]

GROUP BY [S_MFD_ID];

ALTER TABLE [TEMP3] (ADD INDEX [S_MFD_ID_XX] BTREE ([S_MFD_ID]));

UPDATE

(

SELECT

 [1].[MFD_ID],

 [1].[COORD] AS [A],

 [2].[O_GEOM] AS [B]

FROM CALL A_TBL() AS [1]

INNER JOIN [TEMP3] AS [2] ON [1].[MFD_ID] = [2].[S_MFD_ID]

)

SET [A] = [B];

DROP ROOT X;

tjhb

8,958 post(s)
#02-Aug-19 01:26

That's very nice! But somewhat hard to read with no description or commenting. (I predict that you will also find it hard to read in one week or one month.)

I will try to write a simple description of what you have done, before posting my own approach.

Small observation: you have missed a semicolon at the end of the second SELECT INTO statement. However, no error is thrown (I think that is a bug in the parser?), and the following ALTER TABLE statement still works, i.e. adds index [S_XX] (I think it should not).

tjhb

8,958 post(s)
#02-Aug-19 02:37

I will try to write a simple description of what you have done, before posting my own approach.

First a comment. It seems weird to me to use FUNCTION to describe a static TABLE, rather than using VALUE for the purpose. AFAIK they are equivalent, except that I am not sure that the result of the FUNCTION will be cached, while I think the result of VALUE will be. Perhaps I am wrong, and neither can be cached. This could be tested with data.

Anyway, that is very thought-provoking! Your approach makes sense to you and works. Two good things.

Here is my simple description of what your query does.

  • (1) Draw a point for every pixel in the image. (Currently, at the lower left corner. Should it be in the centre? Or the pixel rect? I think the whole pixel rect.) -> TEMP
  • (2) Overlay the pixel points onto containing areas. -> TEMP2
  • (3) For each area, choose the lowest point. -> TEMP3
  • (4) For each area, store the geom for the chosen point from TEMP3.

I think the approach I have in mind will be more efficient. I will describe it before posting the code.

tjhb

8,958 post(s)
#02-Aug-19 02:50

Here is the approach I was thinking of.

  • (1) For each area, draw convex parts.
  • (2) For each convex part, list each tile it touches (by tile rect, keeping tile XY index only).
  • (3) For all distinct listed tiles, list each pixel (pixel XY index and value).
  • (4) For each convex part, choose the lowest/highest touching pixel (by pixel rect).
  • (5) For each area, list the lowest/highest touching pixel (for all of its convex parts).

It is very similar. I was struck by that when I saw your code. It just saves some work. (Though it is very possible that GeomOverlayContainedPar() and related functions already work with convex parts implicitly.)

Dimitri


5,620 post(s)
#02-Aug-19 14:34

Here's another approach. Use the Trace Areas transform (press Edit Query to see the TileTraceAreasPar function it uses) to generate areas populated with the height taken from pixels. You can then intersect those with your areas of interest and then pick out the highest region within each area of interest. May as well let TileTraceAreasPar do the setup work for you.

See the Trace Areas example.

By the way, I tried this with the Mt Hood example DEM that can be downloaded from the Product Downloads page, using a Similarity of 10. It's worth trying just to watch the Preview fill in, after about 10 or 12 seconds, gradually filling in to the summit of Mt. Hood. :-)

The process goes astonishingly quickly, creating about 1.6 million areas. Below is a zoomed in view near the mountain, styled by CB Spectral palette so blue is lower elevations and red is higher. If you wanted to make this more efficient, just flatten the raster so all heights below any in the upper ranges are the same value, to minimize the number of areas created that you don't care about.

Attachments:
mt_hood_trace_areas.png

atrushwo59 post(s)
#02-Aug-19 17:30

Hey Dimitri, Thanks for this alternate solution. Much appreciated.

tjhb

8,958 post(s)
#03-Aug-19 01:26

That is excellent Dimitri. I wouldn't have thought of doing it this way--it would not have been feasible with previous software.

I will remember this and (I'm sure) find new things we can now also do in a similar way.

Horsepower!

atrushwo59 post(s)
#02-Aug-19 15:30

Hey Tim,

Thanks for the detailed response! I agree that my code is undocumented and confusing. I'll try to leave myself some more notes in the future, as yes, I will not understand this once I get my head out of it.

I wasn't aware of VALUES. I'll look into it. I've been using the FUNCTION command almost as a variable to store static things. I'm slowly working on a library of SQL queries to distribute internally. I have an issue where users don't understand all of the places to change table references. Hence me trying to use variables to simplify changes to more complex processes.

I've actually never had a missing semicolon throw an error in M9. I've begun to wonder if they are even necessary. I try to include them, but don't know why I bother sometimes. Thanks for pointing out my missing one.

Perhaps I don't understand how the compiler works, but I'm not sure why convex parts would speed up the process. Would it just be in the scenario where your areas don't touch certain tiles? I'm guessing both processes would be very similar if you had full tile coverage with your areas?

Anyways, I'm interested to hear your thoughts and see your code!

Andrew

tjhb

8,958 post(s)
#03-Aug-19 01:50

Perhaps I don't understand how the compiler works, but I'm not sure why convex parts would speed up the process.

Simpler tests can be used for touching/containing geometry if input geometry is (known to be) convex.

Would it just be in the scenario where your areas don't touch certain tiles?

That is a separate thing. It is the motivation behind my steps 2 and 3. Splitting all pixels for all tiles is likely to be wasteful, if our target areas only touch some of them.

I'm guessing both processes would be very similar if you had full tile coverage with your areas?

If target areas touch all images tiles, or almost all, then yes, prefiltering tiles would not be a benefit (and might be a net cost).

Similarly, if your geometry is already convex, or mainly very simple, then splitting into convex parts would unnecessary at best. Besides, if the built-in functions are good at judging this, and implicitly use convex parts when that helps, but not otherwise, then it would be better to leave it to the engine to decide.

These things are worth testing, with different data scenarios.

tjhb

8,958 post(s)
#10-Aug-19 00:11

Andrew,

I haven't finished working on this.

I seems I have struck a significant bug in the GeomOverlayTouchingPar() function.

I don't yet know if it also applies to the non-Par version, or to related functions, such as GeomOverlayContainedPar() which you used.

[Edit] Same result with non-Par function GeomOverlayTouching(), and GeomOverlayContaining() is also subject to the issue.

I will post about the issue in a separate thread.

In the meantime I would say: check your result against a join, even though a join is slower.

atrushwo59 post(s)
#12-Aug-19 15:49

Hey Tim, Thanks for the heads up. Will check the results and post in your (upcoming) thread if necessary.

tjhb

8,958 post(s)
#14-Aug-19 04:08

Sorry for dragging the chain on this. I have been away. Back to it now (and more to say in this thread of course given the advent of the new functions).

atrushwo59 post(s)
#15-Aug-19 16:37

Hey Tim,

All good, Adam's solution as part of 9.0.169.7 seemed pretty efficient. The new TileGeomToValues exposes the coordinates pretty effectively. It gets around having to break up the surface into pixels and join back. http://www.georeference.org/forum/t148616.17#148628

I checked the GeomOverlayContainedPar() against a join and got the same result. I didn't notice a significant bug, but I very likely could have missed something.

tjhb

8,958 post(s)
#17-Aug-19 06:03

I checked the GeomOverlayContainedPar() against a join and got the same result. I didn't notice a significant bug, but I very likely could have missed something.

I'm very pleased to report that the significant bug was all mine.

I was feeding a test drawing to a geometry overlay, where the test drawing took geometry from the wrong field. Very dumb and very basic, but I managed to overlook it for a long time. I will remember to try avoiding that in future.

It is hard to make the mistake with a spatial join, because of course we write down the name of the fields used in the condition. With a geometry overlay, the name of the fields used are implicit--they are in the drawings' properties... and of course in the code that creates the drawings in the first place (ahem).

Now that TileGeomToValues and its relatives are built in, I don't know whether it's worthwhile posting my query to do the same thing in long form.

It might be. Apart from showing just how nice the new functions are (how much manual coding they can save), the query is an example of manually breaking data into convex parts for (hopefully) an increase in processing speed.

Another reason is that it would be easy-ish to build on the query's infrastructure, to allow writing values from areas to touching pixels.

...But of course, that is likely to be built in before very long as well, which is great.

adamw


8,775 post(s)
online
#17-Aug-19 14:25

Thanks for taking the time to report that the suspected bug was not in our functions, we kind of got worried.

On writing values from geometry to pixels, we will have functions that do this for sure, it's a question of how to package them best. And as frequently happens, there's a different piece of logic this depends on, that is best implemented first. So, the functions will come, but not right away.

tjhb

8,958 post(s)
#17-Aug-19 15:42

I'm sorry it took me so long to be sure--glad to be wrong.

As often is the case, a side-effect of writing a systematic explanation of what I expected to see and what I was seeing, ready for posting as a report, was to bring my own assumptions into focus. When I checked through those assumptions, my error showed up. A well-written bug report often dissolves into nothing.

On the other hand, if I had written a quick report earlier, no doubt you would have spotted my mistake at once. I should have done that.

adamw


8,775 post(s)
online
#12-Aug-19 16:17

Check out 9.0.169.7. It adds the TileGeomToValues function family - these functions take an image and a geom, and return a table of pixel values under the geom. You can then do Min / Max / Median / whatever with the values *and* you have access to pixel coordinates. The geom has to be in the coordinate system of the image.

atrushwo59 post(s)
#12-Aug-19 16:49

Very cool, and very timely Adam.

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