Subscribe to this thread
Home - General / All posts - Profile of a Line - overlay line on raster
GeoKeys12 post(s)
#23-Jun-20 03:11

In Manifold 9 - how do I join a line to a dem? I am ultimately trying to create an elevation profile of a line. I can see this in: http://www.georeference.org/doc/profiles_and_elevations.htm, but not in version 9. Can someone point me in the right direction?

Thanks,

Dimitri


6,364 post(s)
#23-Jun-20 05:56

There's no point and click dialog to do that in 9 yet. You could do that using a mix of dialogs or SQL, but it would be complex. For example, in multiple steps you could transform a line into a sequence of points, then transfer height values from the raster to those points and then visualize the sequence of heights in a charting package.

GeoKeys12 post(s)
#23-Jun-20 08:53

Ok, thanks.

tjhb

9,508 post(s)
#23-Jun-20 09:17

It might be a bit complex, but not very, and moreover it has already been robustly done.

If the OP posts test data, I will post code to do it.

This is not an isolated case, and it is why Manifold should start curating code as well as creating code.

You've made a tool. Keeping everything ego is just silly.

GeoKeys12 post(s)
#24-Jun-20 08:37

What do you need to do this? Just a dem with a polyline in a manifold file?

tjhb

9,508 post(s)
#24-Jun-20 09:08

Yes, and both can be tiny. But as much detail as possible would be good. I.e. a DEM with noticeable relief, and a line that is not very straight (rather than leaving that as a hard case for later).

Sloots

489 post(s)
online
#24-Jun-20 12:25

Here some test data and a query that generates a table of points every 50m along the line with the corresponding height. I used Excel to create the chart. A good start for a better, more comprehensive approach (multiple lines, how to dela with missing values etc.)

Chris

Attachments:
profile.map
profile.png


http://www.mppng.nl/manifold/pointlabeler

GeoKeys12 post(s)
#09-Jul-20 08:04

Hi,

Thanks for your solution. I had some luck with your data, but had a few issues with mine (broken lines, reversed lines, non-sequential branches etc) that I had to deal with. Now I get an answer but I have multiple lines in the same file. I tried to adapt your script to include the line identifiers as:

VALUE @target NVARCHAR = ComponentCoordSystemAutoXY([1475 ELEVATION Raster]);

VALUE @source NVARCHAR = ComponentCoordSystemAutoXY([Lines]);

VALUE @conv TABLE = CALL CoordConverterMake(@target, @source);

SELECT *,

 TileGeomMax([1475 ELEVATION Raster], 0, CoordConvert(@conv, GeomMakePoint([XY]))), GeomMakePoint(XY) as [geom]

INTO

 [points table]

FROM

(

SELECT *,

    SPLIT CALL GeomToCoords(GeomSegmentize([Geom], 50) ) -- Every 50 m a point

FROM 

 [Lines]

THREADS SystemCpuCount()

);

Manifold doesn't like it - Cannot parse query - even though the subquery works, and doesn't have duplicate fields. I am obviously missing something basic. Ideally like to add the distance from the first start of line all the way in 50m increments.

Can you please advise?

GeoKeys12 post(s)
#09-Jul-20 08:49

Got there eventually - works when I provide an alias for the sub query field mfd_id - presume mfd_id is duplicated:

VALUE @target NVARCHAR = ComponentCoordSystemAutoXY([1475 ELEVATION Raster]);

VALUE @source NVARCHAR = ComponentCoordSystemAutoXY([Lines]);

VALUE @conv TABLE = CALL CoordConverterMake(@target, @source);

SELECT xter,

 TileGeomMax([1475 ELEVATION Raster], 0, CoordConvert(@conv, GeomMakePoint([XY]))), GeomMakePoint(XY) as [geom]

INTO

 [points table]

FROM

(

SELECT [Lines].[mfd_id] as xter,

    SPLIT CALL GeomToCoords(GeomSegmentize([Geom], 50) ) -- Every 50 m a point

FROM 

 [Lines]

THREADS SystemCpuCount()

);

adamw


9,470 post(s)
#09-Jul-20 09:02

I think the fix might have been in the switch from * to a specific field, you probably could have left [Lines].[mfd_id] as is without renaming it.

GeoKeys12 post(s)
#09-Jul-20 09:23

I did try that, but I couldn't get it to work without an alias.

Do you have any suggestions on how I determine the distance? In SQL I would do something like: (row_number() over (partition by [XX] Order by [YY]) -1) * 50. The last record wouldn't be accurate, but if the script has nicely divided the lines into 50m interval, it should be fine for the level of accuracy that I am after.

adamw


9,470 post(s)
#10-Jul-20 13:07

You can divide the lines into fixed-length intervals using GeomToCoordsLineSequence. Eg, the inner SELECT with SPLIT could be:

--SQL

...

SELECT [mfd_id]SPLIT CALL GeomToCoordsLineSequence([geom], 0, 1000000, 50)

FROM [lines]

...

The arguments for GeomToCoordsLineSequence are, in order: line, begin distance (I used 0, you can start with, say, 25 meters, the value could be anything), end distance (I used a big number that likely exceeds the length of the line, you can use something else, the value could again be anything, and if it ends up being less than begin, the coordinates would be produced in the reverse order), step (I used 50, the same value as in your query).

Here's what the produced coordinates look like:

See the Decompose to Line Coordinates transform which uses that function.

On MFD_ID, I see why the alias is needed: yes, switching from * to a specific field makes the query parse, but this just fails later on SELECT INTO because the inner SELECT may generate multiple records with the same value of MFD_ID, and MFD_ID has to be unique.

Attachments:
line-coordinates.png

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