georeference.org
Subscribe to this thread
Home - Scripting / All posts - SQL to Create Area "patches" by selecting Adjacent objects in same Drawing
geoworldmaps
12 post(s)
#08-Feb-10 12:51

I have a single drawing of areas. I want to create "patches" where adjacent/touching objects are selected and a separate column indicates a character, or number differentiating each patch from another patch of adjacent/touching areas. I do not want to dissolve/union these adjacent/touching area objects. I want them to remain as separate objects.

Thank you if any one can offer assistance

Eric

tjhb

3,494 post(s)
online
#08-Feb-10 13:14

So you don't want to create patches? You just want to mark patches, right?

geoworldmaps
12 post(s)
#08-Feb-10 13:14

That's correct!

tjhb

3,494 post(s)
online
#08-Feb-10 14:13

Then what's the maximum number of patches you need to allow for?

It doesn't have to be an exact number. A power of ten would be fine.

geoworldmaps
12 post(s)
#08-Feb-10 14:18

More than 100, and possibly up to 10,000..

tjhb

3,494 post(s)
online
#08-Feb-10 15:08

Thanks. Actually it's not necessary to know that. There's a better way.

Name the first query [Islands]. This returns the islands (contiguous patches) in [Drawing].

-- Islands

SELECT [Island]

FROM

    (SELECT UnionAll([Geom (I)]AS [Union]

    FROM [Drawing]

    )

SPLIT BY Islands([Union]AS [Island];

Name the second query [Islands with index]. This calls [Islands] twice, and gives each island an index. The index, for each metric, is the number of metrics less than it. (We don't need to know exactly what "less than" means, provided that metrics in different locations necessarily evaluate to different values, which I believe to be true.)

-- Islands with index

SELECT 

    COUNT([T2].[Island]AS [N],

    [T1].[Island]

FROM

    [Islands] AS [T1]

    LEFT JOIN

    [Islands] AS [T2]

    ON CentroidInner([T1].[Island]) > CentroidInner([T2].[Island])

GROUP BY [T1].[Island];

The third query, [Areas by island], returns each original area, along with the index of the island it belongs to.

-- Areas by island

SELECT

    [D].[Geom (I)],

    [T].[N] AS [Island number]

FROM

    [Drawing] AS [D]

    LEFT JOIN

    [Islands with index] AS [T]

    ON Touches([D].[Geom (I)][T].[Island])

    AND ClipSubtract([D].[Geom (I)][T].[Island]IS NULL;

If you want to put the island index (i.e. patch number) into a column in the original drawing, then the third query can be made into an UPDATE version.

-- Areas by island (UPDATE)

UPDATE

    (SELECT

        [D].[Island][T].[N]

    FROM

        [Drawing] AS [D]

        LEFT JOIN

        [Islands with index] AS [T]

        ON Touches([D].[Geom (I)][T].[Island])

        AND ClipSubtract([D].[Geom (I)][T].[Island]IS NULL

    )

SET [D].[Island] = [N];

If you want to know where the island with each index is, link a drawing from the second query.

geoworldmaps
12 post(s)
#10-Feb-10 09:33

Thank You!

Of course I could have done about 15 steps in ArcGIS, but this is such a great education for me to learn

Spatial SQL in an applied application. Truly appreciate this as it helps me also learn the ins and outs of

querying and geoprocessing in Manifold.

Thanks for the lesson.

artlembo


1,871 post(s)
#10-Feb-10 13:15

Eric,

here is what I think is a similar result, but using VBScript. Perhaps its shorter:

just create an empty drawing called D2

Sub Main

 Set CompS = document.ComponentSet

 Set D = CompS("D")

 Set D2 = CompS("D2")

 D.Copy

 D2.Paste

 Set D_OS = D.ObjectSet

 Set Analz = document.NewAnalyzer

 Analz.Union D2, D2, D2.ObjectSet

 Analz.Decompose D2, D2.ObjectSet

 Set Qry = Document.NewQuery("qry",true)

 Qry.Text = "update D Set island = (select D2id from (select D2.id D2id, D.id Did from D2,D where contains(d2.id,D.id)) where D.id = Did)"

 Qry.Run

End Sub

 

 

Attachments:
islands.map

geoworldmaps
12 post(s)
#10-Feb-10 16:42

It's like dueling banjos out there

artlembo


1,871 post(s)
#10-Feb-10 17:00

ha, ha. Well, Tim is the master when it comes to SQL, and his queries will translate to SQLServer, PostGIS, and Oracle spatial quite well (with a little modification). The code I wrote attempts to leverage the Manifold API which in some cases can be more efficient, but it does in fact lock you in to the Manifold API.

Its great to have both banjos.

0 msec Copyright (C) 2007-2008 Manifold.net. All rights reserved.