georeference.org
Subscribe to this thread
Home - Scripting / All posts - Create rectangles using data from a point
#06-Nov-16 19:45

I am interested in creating rectangles from a table of points.

Each point is the center of a rectangle and contains attributes including width, length, and rotation (direction).

Can anyone get me started on scripting something like this? Or point me to a resource that does something similar (create an object using given coordinates for the center of an object and dimensions)?

Thanks.


cbrown

tjhb
7,048 post(s)
#06-Nov-16 21:52

That is most easily done (and definitely quickest) in Manifold SQL, rather than a script.

Rather than giving hints, better to write an example query, with good comments so you can see how it works.

Can you supply sample data? Even one rectangle is enough (and it can be fake).

#11-Nov-16 13:57

Thanks for giving it some thought.

Attached is a .csv file with some sample data.

Columns in the data are Obj_ID, Width, Distance, Direction, Longitude, Latitude.

Attachments:
Rectangle from point data.csv


cbrown

Sloots

397 post(s)
#07-Nov-16 17:31

This is my version of an SQL query to solve your problem. I have created a table with information about the rectangles (x, y, width, height and angle). Then I wrote two SQL queries, one to show the center points (x, y) just as a reference, and a second one to create the desired rectangles.

Both queries were linked as drawings and shown together in a map. Works great!

Hope this helps.

Cheers,

Chris

Attachments:
rectangles_from_table.map


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

#11-Nov-16 14:09

Thanks! I will take a look at it and let you know how it works for my situation.

Chip


cbrown

#11-Nov-16 21:52

When I use the query on my data (that has lat lon coordinates as opposed to the x-y coordinates in your query) all of the rectangles end up centered at the same point.

Any ideas? Maybe looking at the data I supplied above will help?

Thanks in advance.


cbrown

tjhb
7,048 post(s)
#11-Nov-16 22:36

Columns in the data are Obj_ID, Width, Distance, Direction, Longitude, Latitude.

First you need to say what projection you want to use to measure Width, Distance and Direction.

In other words, what space are these numbers measured in?

If Width and Distance are interpreted in Latitude / Longitude (unprojected space), they will be read as angular units (degrees, minutes or seconds), which I doubt is what you mean. While Direction can be interpreted in Latitude / Longtitude, that will not normally give the same rotation as you would measure on the ground. You probably want Direction to be measured in the same projection as Width and Distance.

From there, Chris's query [qryPoints] needs adjustment to give the points the right projection. Currently they are interpreted as being in Latitude / Longitude by default, which is not what is needed. It's easy to forget this. This line is the problem:

NewPoint([X], [Y]) as geometry,

Also, both queries will need an OPTIONS clause to give the correct projection to the linked drawing(s).

[Added.] It looks as if you're using Width to mean the length of the long side of each rectangle, and Distance to mean the length of the short side. That's fine. But we also need to know how to interpret Direction. Does this mean the rotation of the rectangle (after drawing), or does it mean the direction the rectangle is drawn in (the opposite)?

#14-Nov-16 21:39

Thanks for your thoughts on this.

I left out some details that I thought were unimportant, but share them here.

The data I want to use is yield data from combines that are harvesting corn, soybeans, wheat, etc.

The Distance is the distance traveled by the combine over a given time (usually in 1 second)

The Width is the width of the combine head (usually 12, 20, or 24 feet)

The Direction is the Heading of the combine so it relates to the rotation of the rectangle (0-360 degrees), I think.

We import the data as a table and then copy and paste that table to a points drawing. At that point, the data is unprojected. So it seems that in order for this to work, I need to project the data so the distance and width are interpreted correctly. I wonder how to do that so the measurement can be maintained in feet? Or if that is not possible, I can convert feet to meters.

I think I understand the "options" clause, but could not find a way to identify the proper projection. Perhaps the easiest way is to do that is to use a drawing that is already projected to what I want.

What else should consider?

Again, thanks for your advice and expertise.


cbrown

tjhb
7,048 post(s)
#15-Nov-16 03:07

That's really ideal Chris. Context and purpose always help much more than we realise. But of course we try not to burden people with clutter... we can't really win.

But now you've covered everything, and not too much.

Back soon.

tjhb
7,048 post(s)
#17-Nov-16 23:49

The tricky part here is that we're working in two types of coordinate system at once.

First, the data points (centres) are unprojected (Latitude / Longitude). They're easy to draw of course.

But the Width of each cut and the Distance travelled are in linear units (feet), measured of course from the position of the harvester at that point in time.

That's not so easy. What we ideally need, for each cut, is accurate measurement in a local projection centred on the harvester at that moment. That is, in a series of slightly different local projections, one for each data point. The type of local projection should ideally be equal-area and, since we also care about direction of travel, conformal.

(We could just use a single local projection, for a given field, or for a day's work, or for a county or state... These would be successive good approximations, and probably adequate and sensible for most applications. But I've gone for the ideal.)

This isn't very easy in Manfold 8 SQL (it's easier in upcoming Radian Studio), since we have to fiddle with projections in WKT format (using regular expressions). It's also a bit slow to deal with a different projection for each point.

Here is the query.

Note that it produces output in unprojected space. The rectangles covered by the harvester in real (local) space do not appear quite rectangular when unprojected. That's correct, to be expected. We could equally choose a single local projection for output, and if this is done then real rectangles look like rectangles. (To check: unlink the output drawing, reproject to a local LAEA projection.)

OPTIONS

    CoordSys("Latitude / Longitude"),

    Precision(1e-6)

;

SELECT

    [Obj__Id]-- NB double underscore

    [Width][Distance],

    [Direction],

    [Longitude][Latitude],

    --[baseCS], [localCS],

    [P]-- unprojected

    --[Q], -- in local projection

    Project(

        Rotate(

            ConvertToArea(

                AssignCoordSys(

                    NewLine(

                        NewPoint(CentroidX([Q]) - [dX], CentroidY([Q]) - [dY]),

                        NewPoint(CentroidX([Q]) - [dX], CentroidY([Q]) + [dY]),

                        NewPoint(CentroidX([Q]) + [dX], CentroidY([Q]) + [dY]),

                        NewPoint(CentroidX([Q]) + [dX], CentroidY([Q]) - [dY]),

                        NewPoint(CentroidX([Q]) - [dX], CentroidY([Q]) - [dY])

                        ),

                    WktToCoordSys([localCS])

                    )

                ),

            [Direction]

            ),

        CoordSys("Latitude / Longitude")

        ) AS [Rectangle]

        -- NB rectangular in local LAEA projection,

        -- not exactly rectangular in unprojected space

        -- (normal, to be expected)

FROM

    (SELECT

        [Obj__Id],

        [Width][Distance],

        [dX][dY],

        [Direction],

        [Longitude][Latitude],

        [baseCS][localCS],

        [P]-- unprojected

        Project([P], WktToCoordSys([localCS])) AS [Q] -- in local LAEA projection

    FROM

        (SELECT

            [Obj__Id]-- NB double underscore

            [Width][Distance]-- in feet

            [Width] * 12 * 25.4 / 1e3 / 2 AS [dX]-- x offset in metres

            [Distance] * 12 * 25.4 / 1e3 / 2 AS [dY]-- y offset in metres

            [Direction],

                -- bearing, assumed relative to true north

            [Longitude][Latitude],

                -- centre of rectangle harvested in unit time

                -- (question: should this instead be the centre

                -- of the upper edge of the harvested rectangle?)

            [baseCS],

            RegExp(

                RegExp(

                    [baseCS],

                    [getLong][newLong] & CStr([Longitude])

                    ),

                [getLat][newLat] & CStr([Latitude])

                ) AS [localCS],

            NewPointLatLon([Longitude][Latitude])

--            AssignCoordSys(

--                NewPoint([Longitude], [Latitude]),

--                CoordSys("Latitude / Longitude")

--                ) -- equivalent

                AS [P] -- unprojected

        FROM

            [Rectangle from point data]

            CROSS JOIN

            (VALUES (

                -- default LAEA projection centred at (0, 0) as WKT

                CoordSysToWkt(CoordSys("Lambert Azimuthal Equal Area")),

                -- search strings for longitude, latitude parameters in WKT

                "\x22Longitude_Of_Center\x22,(\d*\.*)?\d+",

                "\x22Latitude_Of_Center\x22,(\d*\.*)?\d+"

                -- replacement strings (without number)

                Chr(34) & "Longitude_Of_Center" & Chr(34) & ","

                Chr(34) & "Latitude_Of_Center" & Chr(34) & ","

                )

            NAMES ([baseCS][getLong][getLat][newLong][newLat])

            )

        )

    )

;

There's a question in-line. I wondered whether the rectangle for each cut should be drawn with the data point in their centre, or with the data point at their furthest edge. It depends when and where the harvest for each cut is measured. I've used centres, since that's what was mentioned.

Attachments:
20161118 Draw unprojected rectangles measured in local LAEA projections.map
Draw rectangles.txt

tjhb
7,048 post(s)
#18-Nov-16 01:26

It only occurred to me afterwards that, because each point Q is the same as the origin of its local projection, we don't need it. It's always (0, 0).

So we can get rid of point Q, and draw each rectangle using just the offsets dX and dY.

This simpler version is also faster to link.

OPTIONS

    CoordSys("Latitude / Longitude"),

    Precision(1e-6)

;

SELECT

    [Obj__Id]-- NB double underscore

    [Width][Distance],

    [Direction],

    [Longitude][Latitude],

    --[baseCS], [localCS],

    [P]-- unprojected

    Project(

        Rotate(

            ConvertToArea(

                AssignCoordSys(

                    NewLine(

                        NewPoint(-[dX], -[dY]),

                        NewPoint(-[dX], +[dY]),

                        NewPoint(+[dX], +[dY]),

                        NewPoint(+[dX], -[dY]),

                        NewPoint(-[dX], -[dY])

                        ),

                    WktToCoordSys([localCS])

                    )

                ),

            [Direction]

            ),

        CoordSys("Latitude / Longitude")

        ) AS [Rectangle]

        -- NB rectangular in local LAEA projection,

        -- not exactly rectangular in unprojected space

        -- (normal, to be expected)

FROM

    (SELECT

        [Obj__Id],

        [Width][Distance],

        [dX][dY],

        [Direction],

        [Longitude][Latitude],

        [baseCS][localCS],

        [P]

    FROM

        (SELECT

            [Obj__Id]-- NB double underscore

            [Width][Distance]-- in feet

            [Width] * 12 * 25.4 / 1e3 / 2 AS [dX]-- x offset in metres

            [Distance] * 12 * 25.4 / 1e3 / 2 AS [dY]-- y offset in metres

            [Direction],

                -- bearing, assumed relative to true north

            [Longitude][Latitude],

                -- centre of rectangle harvested in unit time

                -- (question: should this instead be the centre

                -- of the upper edge of the harvested rectangle?)

            [baseCS],

            RegExp(

                RegExp(

                    [baseCS],

                    [getLong][newLong] & CStr([Longitude])

                    ),

                [getLat][newLat] & CStr([Latitude])

                ) AS [localCS],

            NewPointLatLon([Longitude][Latitude])

--            AssignCoordSys(

--                NewPoint([Longitude], [Latitude]),

--                CoordSys("Latitude / Longitude")

--                ) -- equivalent

                AS [P] -- unprojected

        FROM

            [Rectangle from point data]

            CROSS JOIN

            (VALUES (

                -- default LAEA projection centred at (0, 0) as WKT

                CoordSysToWkt(CoordSys("Lambert Azimuthal Equal Area")),

                -- search strings for longitude, latitude parameters in WKT

                "\x22Longitude_Of_Center\x22,(\d*\.*)?\d+",

                "\x22Latitude_Of_Center\x22,(\d*\.*)?\d+"

                -- replacement strings (without number)

                Chr(34) & "Longitude_Of_Center" & Chr(34) & ","

                Chr(34) & "Latitude_Of_Center" & Chr(34) & ","

                )

            NAMES ([baseCS][getLong][getLat][newLong][newLat])

            )

        )

    )

;

Attachments:
Draw rectangles b.txt
Image.png

#18-Nov-16 18:51

Wow! That is going to take me a while to digest.

Thanks. I am eager to try it and the image looks like it works well!

As far as your question about whether the point is the center or the upper edge of the rectangle: That is a good question and I will need to do a little investigation to figure that out. I have always assumed that it is the center point, but once you mentioned it, it also makes sense that it would be the upper edge of the rectangle.

I will give it a go and report back!

-Chip


cbrown

tjhb
7,048 post(s)
#18-Nov-16 19:54

Thanks for thinking about that question (I didn't manage to put it very clearly).

I think your intuition is right (mine too), that the upper edge (leading edge) of the rectangle should be used. In other words, that the rectangle should be drawn behind its report point, not centred on it.

Apart from inituition, the reason I think this looks right is that if the rectangles are drawn this way, then there are fewer apparent gaps and overlaps in coverage, that is, areas which were apparently not covered, or covered twice, by the harvester. (There are almost none now.)

To see that, contrast attached images "Map image geographic b.png" (using centres) with "Map image geographic c.png" (using the centre of each leading edge).

Likewise, "Map image LAEA b.png" versus "Map image geographic c.png". These two are made from a map using a Lambert Azimuthal Equal Area projection centred at -77.5, 36.5 (roughly the centre of this data). In this (average) local projection you can also see that the rectangles appear as rectangles.

Here is the query adjusted to use the centre of the leading edge. (Differences in headnote.)

-- Generally follows version b

-- but assumes that the report point for each cut

-- follows the cut

-- So the harvested rectangle is drawn behind the report point,

-- rather than centred on it as for versions a and b

-- See different treatment of [dY]

-- Requires RotateAbout() rather than Rotate()

-- (Also a cosmetic change: AssignCoordSys() is applied to area rectangle

-- rather than the line)

OPTIONS

    CoordSys("Latitude / Longitude"),

    Precision(1e-6)

;

SELECT

    [Obj__Id]-- NB double underscore

    [Width][Distance],

    [Direction],

    [Longitude][Latitude],

    --[baseCS], [localCS],

    [P]-- unprojected

    Project(

        RotateAbout(

            AssignCoordSys(

                ConvertToArea(

                    NewLine(

                        NewPoint(-[dX], -[dY]),

                        NewPoint(-[dX], 0),

                        NewPoint(+[dX], 0),

                        NewPoint(+[dX], -[dY]),

                        NewPoint(-[dX], -[dY])

                        )

                    ),

                WktToCoordSys([localCS])

                ),

            AssignCoordSys(

                NewPoint(0, 0), 

                WktToCoordSys([localCS])

                ),

            [Direction]

            ),

        CoordSys("Latitude / Longitude")

        ) AS [Rectangle]

        -- NB rectangular in local LAEA projection,

        -- not exactly rectangular in unprojected space

        -- (normal, to be expected)

FROM

    (SELECT

        [Obj__Id],

        [Width][Distance],

        [dX][dY],

        [Direction],

        [Longitude][Latitude],

        [baseCS][localCS],

        [P]

    FROM

        (SELECT

            [Obj__Id]-- NB double underscore

            [Width][Distance]-- in feet

            [Width] * 12 * 25.4 / 1e3 / 2 AS [dX],

                -- lateral offset in metres

            [Distance] * 12 * 25.4 / 1e3 AS [dY],

                -- distance travelled in metres

                -- (assumes report follows harvest)

            [Direction],

                -- bearing, assumed relative to true north

            [Longitude][Latitude],

                -- report point for rectangle harvested in unit time

            [baseCS],

            RegExp(

                RegExp(

                    [baseCS],

                    [getLong][newLong] & CStr([Longitude])

                    ),

                [getLat][newLat] & CStr([Latitude])

                ) AS [localCS],

            NewPointLatLon([Longitude][Latitude])

--            AssignCoordSys(

--                NewPoint([Longitude], [Latitude]),

--                CoordSys("Latitude / Longitude")

--                ) -- equivalent

                AS [P] -- unprojected

        FROM

            [Rectangle from point data]

            CROSS JOIN

            (VALUES (

                -- default LAEA projection centred at (0, 0) as WKT

                CoordSysToWkt(CoordSys("Lambert Azimuthal Equal Area")),

                -- search strings for longitude, latitude parameters in WKT

                "\x22Longitude_Of_Center\x22,(\d*\.*)?\d+",

                "\x22Latitude_Of_Center\x22,(\d*\.*)?\d+"

                -- replacement strings (without number)

                Chr(34) & "Longitude_Of_Center" & Chr(34) & ","

                Chr(34) & "Latitude_Of_Center" & Chr(34) & ","

                )

            NAMES ([baseCS][getLong][getLat][newLong][newLat])

            )

        )

    )

;

Revised .map file and query also attached.

Attachments:
Draw rectangles c.txt
Map image geographic b.png
Map image geographic c.png
Map image LAEA b.png
Map image LAEA c.png

tjhb
7,048 post(s)
#18-Nov-16 20:00

Revised .map file.

Attachments:
20161118 Draw unprojected rectangles measured in local LAEA projections.map

#30-Nov-16 01:01

Again, many thanks for working on this. This solution is certainly not one I would have ever figured out.

Harvest season is over and now we have lots of data to process and I think this will enhance what we are putting out there.


cbrown

tjhb
7,048 post(s)
#30-Nov-16 02:09

Isn't it interesting, when a question that can be put so simply in human terms should have such a tricky answer. I love this sort of problem.

Thanks for wading through the detail.

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