Subscribe to this thread
Home - General / All posts - Dot Density map
apo108 post(s)
#04-Jul-20 21:45

Based on the discussion of an older Thread (http://www.georeference.org/forum/t136595.9) I just composed a small query to generate dot density map

The approach is based on the triangulation of the complex geometries cutting the population into parts based on the area of the triangles and then generating the random distribution of the points in the triangles.

There is a side-effect related to the number of triangles needed to compose complex forms generating a lot of those oversplitting the number of points by triangle and depending on the rounding method used (floor, round, ceil) is over or under sampling the points. Quite efficient even with millions of points. I think the code is self explicit, just because the name of the function in M9 are on my point of view explicit. For sure this query can be improved, just a first test

The commented lines are used in case of dynamic drawing but with 10mio points I prefered an INTO

The popfactor variable is used to generalized the number of person by point, in this example 1

code

-- $manifold$

VALUE @popfactor INT32 = 1;

--TABLE CALL TableCacheIndexGeoms(

--  (

SELECT     [mfd_id]/1 as mfdid,

TriID,

GeomMakePoint(

    VectorMakeX2(

    VectorValue(A,0)+[value]*(VectorValue(B,0)-VectorValue(A,0))+((1-[value])*[value 2])*(VectorValue(C,0)-VectorValue(A,0))

    ,

    VectorValue(A,1)+[value]*(VectorValue(B,1)-VectorValue(A,1))+((1-[value])*[value 2])*(VectorValue(C,1)-VectorValue(A,1))

    ) 

) as [Geom]

INTO PtDensity

FROM (

    select [mfd_id],

    TriID,

    A,

    B,

    C,

    value,

    SPLIT CALL ValueSequenceRandom(1, value)

    FROM ( 

        select [mfd_id],

        TriID,

        A,

        B,

        C,

        SPLIT CALL ValueSequenceRandom(TriPtsNbr, TriID)

 

        FROM (

            select [mfd_id],

            TriID,

            avg([TriPtsNbr]) as TriPtsNbr,

            max(case when [coord] = 0 then [xy] else NULL end) as A, 

            max(case when [coord] = 1 then [xy] else NULL end) as B,

            max(case when [coord] = 2 then [xy] else NULL end) as C

            FROM (

 

                select [mfd_id], 

                [branch] as TriID,

                [pop],

                GeomArea([Value],0)/10000 as Triarea,

                GeomArea([Value],0)/10000/[Area] as Tripart,

                round([pop]*GeomArea([Value],0)/10000/[Area]/@popfactor) as TriPtsNbr,

                [value] as [TriGeom],

                SPLIT CALL GeomToCoords([VALUE])

                from (

 

                    SELECT [mfd_id], 

                    [pop],

                    geomarea([Geom],0)/10000 as [Area],

                    [pop] / (geomarea([Geom],0)/10000) as Density,

                    split call GeomtoBranches(GeomTriangulate([Geom],0))

                    FROM [Drawing]

 

                )

            )

            group by [mfd_id],TriID

        )

    )

)

--)

--  , True

--)

;

apo108 post(s)
#04-Jul-20 22:10

The rendering is efficient as show on the fist joined example and might also integrated thematic colouring.

Zooming we can see a side effect considering always the same point A from the triangle. the high densities tends to concentrate the points to the A corner of the triangle. I trick against that would be to better manage the random or randomly select the A point among the 3 corners of the triangle

Attachments:
dotdensmap.png
dotdensmap_sideeffect.png

tjhb

9,320 post(s)
#04-Jul-20 23:01

Very nice!

If you like you can remove the step to

SPLIT CALL GeomToCoords([VALUE])

and the subsequent grouping and pivot. It would be enough just to address the triangle vertices by index.

apo108 post(s)
#04-Jul-20 23:09

Now you say it it seems obvious, thanks I'll try it.

a.

apo108 post(s)
#07-Jul-20 09:53

the new version based on the inputs suggested by Tim

code

-- $manifold$

VALUE @popfactor INT32 = 1;

--TABLE CALL TableCacheIndexGeoms(

--  (

SELECT     [mfd_id]/1 as mfdid,

TriID,

GeomMakePoint(

    VectorMakeX2(

    VectorValue(A,0)+[value]*(VectorValue(B,0)-VectorValue(A,0))+((1-[value])*[value 2])*(VectorValue(C,0)-VectorValue(A,0))

    ,

    VectorValue(A,1)+[value]*(VectorValue(B,1)-VectorValue(A,1))+((1-[value])*[value 2])*(VectorValue(C,1)-VectorValue(A,1))

    ) 

) as [Geom]

INTO PtDensitypopjob_

FROM (

    select [mfd_id],

    TriID,

    A,

    B,

    C,

    value,

    SPLIT CALL ValueSequenceRandom(1, value)

    FROM ( 

        select [mfd_id],

        TriID,

        A,

        B,

        C,

        SPLIT CALL ValueSequenceRandom(TriPtsNbr, TriID)

 

        FROM (

 

                select [mfd_id], 

                [branch] as TriID,

                [pop],

                GeomArea([Value],0)/10000 as Triarea,

                GeomArea([Value],0)/10000/[Area] as Tripart,

                round([pop]*GeomArea([Value],0)/10000/[Area]/@popfactor) as TriPtsNbr,

                [value] as [TriGeom],

                GeomCoordXY([VALUE], 0) as A, 

                GeomCoordXY([VALUE], 1) as B,

                GeomCoordXY([VALUE], 2) as C

                from (

 

                    SELECT [mfd_id], 

                    [pop18]+[j] as [pop],

                    geomarea([Geom],0)/10000 as [Area],

                    split call GeomtoBranches(GeomTriangulate([Geom],0))

                    FROM [Drawing]

 

                )

        )

    )

)

--)

--  , True

--)

;

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