Subscribe to this thread
Home - General / All posts - Extracting lat lon from Geom (I)
ranger.sean109 post(s)
#11-Feb-11 09:01

Okay so I have a drawing which has a number of areas. What I somehow need to do is extract the lat lon for each vertex of each area (some have dozens). I can recall a thread that touched on this but can't find it. Can anyone help?

tjhb
10,094 post(s)
#11-Feb-11 09:11

I don't know what your underlying objective is but something like this.

OPTIONS COORDSYS("Drawing" AS COMPONENT);

SELECT

    [ID][Point]

    CentroidX([Point]AS [Longitude]

    CentroidY([Point]AS [Latitude]

FROM

    (SELECT

        [ID]

        Project([Vertex]COORDSYS("Latitude / Longitude")) AS [Point]

    FROM [Drawing]

    WHERE IsArea([ID])

    SPLIT BY Coords([Geom (I)]AS [Vertex]

    )

;

ranger.sean109 post(s)
#11-Feb-11 10:21

Tim that is awesome! Worked a treat.

One more question and I guess it's more of a general SQL question: I now have a table with various distinct lat lon pairs for each ID. Some IDs are repeated over literally dozens of rows according to how complex the geometry of each area is. How would I go about concatenating the data for each ID so that there is only one row for each ID. For example I'm picturing a comma separated series of values in each column. Doable?

mlinth
447 post(s)
#11-Feb-11 10:35

there's a script that does something similar here.

tjhb
10,094 post(s)
#11-Feb-11 10:57

Or you could try this instead of using SPLIT BY. (I think this is probably less like what you want.)

SELECT

    [ID],

    CAST(

        CGeomWKB(

            Project(

                [Geom (I)],

                COORDSYS("Latitude / Longitude")

                )

            )

        AS TEXT

        ) AS [WKT]

FROM [Drawing];

ranger.sean109 post(s)
#11-Feb-11 18:14

Thanks guys.

At face value the script cited by Martin in http://www.georeference.org/forum/t104463.3 would seem to do exactly what I need but I can't get it to work.

tosborn

231 post(s)
#11-Feb-11 22:34

Here is a vbscript scripting approach to 1) get the vertices and 2) format them as you wanted. If you want the coordinates to be in Lat-Lon, the drawing needs to be in that "projection" prior to running the script. I put the results in a comment component call Vertex Report.

Sub Main

    set report = document.newcomments("Vertex Report")

    set thedrawing = document.componentset("Drawing")

    set theobjset = thedrawing.objectset

    for each myarea in theobjset

        report.addtext myarea.id & ": "

        set thebrnchset = myarea.geom.branchset

            i=0

            for each brnch in thebrnchset

                report.addtext "Branch " & i & ": "

                set thepointset = myarea.geom.branchset(i).pointset

                for each vertex in thepointset

                    report.addtext "(" & vertex.x & ", " & vertex.y & ") "

                next

            i=i+1

            next

        report.addtext vbcrlf

    next

End Sub

mlinth
447 post(s)
#11-Feb-11 22:55

Nice :) - and shows that Manifold offers many ways to the same goal!

FWIW here's a version of the script, project with components attached. A good oppportunity to test the relative speeds of SQL and scripting...

M

Sub Main

' This script displays a list of vertices with the lat / lon pairs in one column

' first get all the vertices in a "long list"

' Then loop through the list, and summarise the output into one record per id

' The project needs a query called "query", a query called "actionQuery" 

' and a table called "Table" with columns ID (integer) and coordinates (text)

' And, of course, the drawing to split

 Set doc = Application.ActiveDocument ' Get the doc

 Set cs = doc.ComponentSet ' Get the Component Set

' empty the output table

cs("Table").clear

 

set theQuery = cs("Query") ' the split query

set actionQuery = cs("ActionQuery") ' this is the insert query for each id

' First run the "split" query to get the vertices

' This is Tim's query, I added sorting by id and longitude

theText = "OPTIONS COORDSYS(""Drawing"" AS COMPONENT);"

theText = theText & "SELECT "

theText = theText & "[ID], [Point], "

theText = theText & "    CentroidX([Point]) AS [Longitude], "

theText = theText & "    CentroidY([Point]) AS [Latitude] "

theText = theText & "FROM "

theText = theText & "    (SELECT "

theText = theText & "        [ID], "

theText = theText & "        Project([Vertex], COORDSYS(""Latitude / Longitude"")) AS [Point] "

theText = theText & "    FROM [Drawing]"

theText = theText & "    WHERE IsArea([ID]) "

theText = theText & "    SPLIT BY Coords([Geom (I)]) AS [Vertex] "

theText = theText & "    ) "

theText = theText & "ORDER BY "

theText = theText & " [ID],Centroidx(Point);"

theQuery.text = theText

theQuery.RunEx TRUE ' run the query

Set theTable = theQuery.Table ' the querie's table

 Set theRecords = theTable.RecordSet '... and the recordset

recordCount = 0

' track the previous id. As we loop through the records, write a record when this changes

oldID = "" ' the id from the previous record

theValue = "" ' A list of all the lon / lat pairs

newID = "" ' the ID from the current record

newValue = "" ' the lon / lat pair from the current record

FOR EACH r in theRecords  ' loop through the records

 recordCount = recordCount + 1 ' increment the record counter

 

 newID = r.datatext("ID") ' get the new ID, treat it as text in the script

  newValue = r.dataText("Longitude") & " " & r.dataTExt("Latitude") ' concatenate lon and lat

 

  IF newID <> oldID THEN ' The ID has changed

  ' IF it is not the first record, write the old record

  IF recordcount > 1 THEN

  actionQuery.text = "insert into [Table] values(" & oldID & " ,""" & theValue & """)" 

   actionQuery.RunEx TRUE

  END IF

 

  ' reset the variables

  oldID = newID

  theValue = newValue

  ELSE ' append the next value

  theValue = theValue & ", " & newValue 

  END IF

 

 NEXT

' write the last record

 actionQuery.text = "insert into [Table] values(""" & oldID & """ ,""" & theValue & """)" 

 actionQuery.RunEx TRUE

' display the table

 cs("Table").open

END SUB

Attachments:
vertices.map

Graeme

990 post(s)
#07-Jul-16 02:16

I had someone seek the same thing so found this thread - tjhb's first query works great, thanks.

There's an additional requirement though and I'm struggling. For each [ID], is it possible to add a column with an incrementing integer value from 1 to n for each coordinate? Can't seem to find an operator equivalent to arithmetic series in transform toolbar.

tjhb
10,094 post(s)
#07-Jul-16 03:26

There are two ways, both with minor drawbacks.

One is to split the original line (or boundary) into two parts at each coord, then take the first part and count its coords. The downside here is that IntersectLine also resolves any self-intersections in the line [and also removes any redundant coords, messing up the index]. If there aren't any then that's not a problem. Otherwise a lengthy workaround is possible.

The second way is simply to join each metric to a simple table of values, 0 to N, where N < the number of coords for the metric, and use these values to look up each coord. The downside here is having to know in advance a reasonable maximum number of coords for all metrics in the set. Usually not tricky, so if that's no bother then this method is easier.

Graeme

990 post(s)
#07-Jul-16 03:43

The data should be clean. I've attached a small sample including your earlier query.

The full data set would have some branched areas and I'm not sure whether they'd pose an issue - for example, casting as WKT for some destination environments requires updating branched records to "multipolygon" with an additional pair of brackets.

There's a handy sample map with good scripts to add a sequential integer, but beyond me to add the code to restart from 1 for each ID.

Attachments:
Rows of area coords with sequential number for each test.map

tjhb
10,094 post(s)
#07-Jul-16 23:21

Graeme,

The full data set would have some branched areas and I'm not sure whether they'd pose an issue - for example, casting as WKT for some destination environments requires updating branched records to "multipolygon" with an additional pair of brackets

I can't quite work this out, sorry for being thick. What format do you need?

Breaking everything into individual points with an index will "dissolve" branches, listed essentially in arbitrary order. But we could break metrics into <object-index>:<branch-index>:<point> or ...<x>,<y> if you like.

But do you ultimately want instead want WKT? That's what I'm confused about there I suppose.

Looking at the sample data, the number of coords is uniformly small (< 200), so there's no problem using a simple join. The same can be done for branches, up to some approximate maximum. (But subject to the above, i.e. whether you want WKT instead.)

Graeme

990 post(s)
#08-Jul-16 01:34

Sorry for the confusion Tim, I used the WKT example as an illustration of the need to explicitly accommodate branched areas if generating WKT destined for certain end user applications - outside Manifold. So no, WKT is not the objective here.

The full data will have instances where there are many more than 200 coords. It is postcode areas, the example covers a small part of metro Melbourne, Australia, so postcodes tend to be small geographically. Moving to remote areas especially coastal ones with lots of islands, the coord number might be up to 50000 with over 1500 branches. In a pragmatic sense, most users of postcode and similar demographic geography, are focused on more densely populated urban areas, leaving the option of filtering out large areas, if required.

I can't say whether

Breaking everything into individual points with an index will "dissolve" branches

that would be problematic for the potential end user - I still don't have details, just a simple spreadsheet of what table structure/format they (believe they) require. A case (not that uncommon) of someone suddenly thinking geography is potentially useful to them and hoping to just push "it" into a system designed for something different. But that approach may well work just fine. I wanted to cover the bases by being able to explicitly address each branch, if required. Sorry it's still a bit vague..

mdsumner


4,260 post(s)
#08-Jul-16 11:12

I want this too, can probably do it just haven't had a closer look again since last discussed!


https://github.com/mdsumner

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