ST_TriangularGrid and ST_HexagonalGrid not working in spatial

[Shane Woodard]

I was wondering if anyone has had trouble creating triangle or hexagonal grids with spatialite gui version 1.8.0-devel?  I can use ST_SquareGrid fine, but either triangle or hexagonal both return null values...

Any suggestions would be greatly appreciated.  I have attached some screenshots:

the spatialite-gui version:

 the sql for the functional square and the results:

the sql I am trying to use for the hexagon which returns a null

Thank you again for any help or pointers!!!!

Shane

[Jukka Rahkonen]

This works for me:

SELECT
ST_TriangularGrid(ST_GeomFromText('POLYGON (( 0 0, 0 120, 160 130, 140 0, 0 0 ))'),25,0)

As well as this:
SELECT
ST_HexagonalGrid(ST_GeomFromText('POLYGON (( 0 0, 0 120, 160 130, 140 0, 0 0 ))'),25,0)

Try these. If they work the issue may have something to do with the 'Smith' geometry. If it ain't very large paste it as WKT for testing.

-Jukka Rahkonen-

[Shane Woodard]

Jakka,

Thank you for the response, and your version does work in my setup as well.  The Smith shape is rather large (6258 vertices), so I zipped it up and attached it if you would like to look at it.  I'll do some more work on my side as well know that I know it does work.

Thank you for your help!!! 

[Shane Woodard]

also as wkt attachment.

[Shane Woodard]

Well, still no luck.  I tried simplifying it down to Jakka's example and used the following two commands.  

I used envelope to create a general polygon (of the Smith county shape I was originally trying to use) and then converted it to WKT.  

I did this as both SRID 4269 (the first example) and then converted it to SRID 3857 as I was wondering if there was some issue with lat/long.  Neither of the below examples work for hexagonal or triangular grid, but both work for square.....

SELECT 

ST_HexagonalGrid(ST_GeomFromText('

POLYGON((

-95.594541 32.135399, 

-94.985265 32.135399,

-94.985265 32.687026, 

-95.594541 32.687026,

-95.594541 32.135399))

'), .3) as geom

 

 

SELECT 

ST_TriangularGrid(ST_GeomFromText(' 

POLYGON((

-10641535 3781097,

-10573711 3781097,

-10573711 3853835,

-10641535 3853835,

-10641535 3781097))

'), 10000) as geom

On Tuesday, June 2, 2015 at 12:06:43 PM UTC-5, Jukka Rahkonen wrote:

[Shane Woodard]

It appears to be an issue with negative numbers....  I can take the statements below and use only positive numbers and everything works as expected...

[Jukka Rahkonen]

I made a similar test with

SELECT
st_hexagonalgrid(st_geomfromtext('POLYGON (( -95.60000642169064 32.69515753782014, -94.98823776665468 32.60149364718705, -94.97220520879856 32.12895509984892, -95.48355942252516 32.12811128101439, -95.61013224770502 32.48420282918705, -95.60000642169064 32.69515753782014 ))'),0.3);

No luck for me either. I would make one more test by turning the negative longitudes into positive but anyway this feels like a bug somewhere.

-Jukka-

[Brad Hards]

Interestingly it does seem to work for the square grid, which might help with
debugging this.

Brad

[Brad Hards]

On Tue, 2 Jun 2015 11:59:50 AM Shane Woodard wrote:
> It appears to be an issue with negative numbers.... I can take the
> statements below and use only positive numbers and everything works as
> expected...

Looks like the problem is based on an assumed origin of (0,0). As a
workaround, if you explicitly set the origin to the bottom left point, then it
should give you a solution.

SELECT AsText(st_triangulargrid(st_geomfromtext('POLYGON (( -95.60000642169064

32.69515753782014, -94.98823776665468 32.60149364718705, -94.97220520879856
32.12895509984892, -95.48355942252516 32.12811128101439, -95.61013224770502

32.48420282918705, -95.60000642169064 32.69515753782014 ))'),0.3, 0,
MakePoint(-95.61013224770502, 32.12895509984892)));
MULTIPOLYGON(((-95.610132 32.128955, -95.310132 32.128955, -95.460132
32.388763, -95.610132 32.128955)), ((-95.460132 32.388763, -95.310132
32.128955, -95.160132 32.388763, -95.460132 32.388763)), ((-95.310132
32.128955, -95.010132 32.128955, -95.160132 32.388763, -95.310132 32.128955)),
((-95.160132 32.388763, -95.010132 32.128955, -94.860132 32.388763, -95.160132
32.388763)), ((-95.010132 32.128955, -94.710132 32.128955, -94.860132
32.388763, -95.010132 32.128955)), ((-95.760132 32.388763, -95.460132
32.388763, -95.610132 32.64857, -95.760132 32.388763)), ((-95.610132 32.64857,
-95.460132 32.388763, -95.310132 32.64857, -95.610132 32.64857)), ((-95.460132
32.388763, -95.160132 32.388763, -95.310132 32.64857, -95.460132 32.388763)),
((-95.310132 32.64857, -95.160132 32.388763, -95.010132 32.64857, -95.310132
32.64857)), ((-95.160132 32.388763, -94.860132 32.388763, -95.010132 32.64857,
-95.160132 32.388763)), ((-95.610132 32.64857, -95.310132 32.64857, -95.460132
32.908378, -95.610132 32.64857)), ((-95.460132 32.908378, -95.310132 32.64857,
-95.160132 32.908378, -95.460132 32.908378)), ((-95.310132 32.64857,
-95.010132 32.64857, -95.160132 32.908378, -95.310132 32.64857)))

Brad

[Brad Hards]

On Wed, 3 Jun 2015 11:28:05 AM Brad Hards wrote:
> On Tue, 2 Jun 2015 11:59:50 AM Shane Woodard wrote:
> > It appears to be an issue with negative numbers.... I can take the
> > statements below and use only positive numbers and everything works as
> > expected...
>
> Looks like the problem is based on an assumed origin of (0,0). As a
> workaround, if you explicitly set the origin to the bottom left point, then
> it should give you a solution.

I spoke too soon. It looks like there is at least one more bug, which is
actually shown by one of the test cases:

spatialite> SELECT AsText(TriangularGrid(GeomFromText('POLYGON((0.2 0.2, 2.2
0.2, 2.2 1.2, 0.2 1.2, 0.2 0.2))', 4326), 1.0));
MULTIPOLYGON(((0 0.866025, 1 0.866025, 0.5 1.732051, 0 0.866025)), ((0.5
1.732051, 1 0.866025, 1.5 1.732051, 0.5 1.732051)), ((1 0.866025, 2 0.866025,
1.5 1.732051, 1 0.866025)), ((1.5 1.732051, 2 0.866025, 2.5 1.732051, 1.5
1.732051)), ((2 0.866025, 3 0.866025, 2.5 1.732051, 2 0.866025)))
spatialite> SELECT ST_MinY(TriangularGrid(GeomFromText('POLYGON((0.2 0.2, 2.2
0.2, 2.2 1.2, 0.2 1.2, 0.2 0.2))', 4326), 1.0));
0.866025403784439

That can't be right - we are supposed to be covering the area (exactly), and
the Y extent doesn't go far enough down.

Brad

[Shane Woodard]

FYI for the next person to run into this.... :-)

Here is what I can up with for a final solution at this time:

select transform(ST_HexagonalGrid(transform(buffer(geometry,.01), 3081), 1000),4269) as geometry

from tx_county

where name = "Smith";

Since the negative numbers seem to cause the issue, I need to reproject into a SRID that covers my area without negative numbers.  In this case I used SRID 3081 which is the NAD83 / Texas State Mapping System.  I noticed when I did that however, there was a small sliver of the south side of the county that was not covered.  I assume that is a reprojection issue and in this case, as long as I got a grid to cover it, I can work with it.  So, working from inside the selection out:

1. buffer the county geometry (in SRID 4269) by .01 to ensure that the reprojection covers the entire county
   
2. transform the buffered geometry into SRID 3081
3. create the grid (which will be in SRID 3081)
4. transform the grid back to 4269 (the original shapes SRID)

This yields a WKT of the final geometry that looks like this (changed the buffer in the statement above to 30000 to keep the WKT smaller):

MULTIPOLYGON(((-95.325687 32.141601,

 -95.178556 31.901227, -94.861101 31.889027,

 -94.689201 32.117142, -94.835417 32.358046,

 -95.154454 32.370305, -95.325687 32.141601)),

 ((-95.792765 32.392478, -95.644049 32.152662,

 -95.325687 32.141601, -95.154454 32.370305,

 -95.302261 32.610651, -95.622215 32.621764,

 -95.792765 32.392478)), ((-95.302261 32.610651,

 -95.154454 32.370305, -94.835417 32.358046,

 -94.662594 32.586076, -94.809475 32.826947,

 -95.130111 32.839263, -95.302261 32.610651)),

 ((-95.771621 32.861542, -95.622215 32.621764,

 -95.302261 32.610651, -95.130111 32.839263,

 -95.278598 33.079566, -95.60016 33.090732,

 -95.771621 32.861542)))

Finally, assuming I save these results to a new table, I can then use the "separate elementary geometries" command in the gui (by right clicking on the geometry column in the new table) , or the .elemgeo macro in the command line, to break out the individual polygon's for my next step in processing data.

Thank you to Jakka who gave the the pointers to help me get to the this solution!

[Brad Hards]

On Wed, 3 Jun 2015 06:49:22 AM Shane Woodard wrote:
> I noticed
> when I did that however, there was a small sliver of the south side of the
> county that was not covered.

That part appears to be the same bug I saw - is it of the order a few hundred
metres (maximum size: of just less than the grid spacing)?

We'll come up with a proper fix, but it isn't a really simple problem...

Brad

[Shane Woodard]

Brad,

I ran the grid to be 5000 meters with the following command:

select transform(ST_HexagonalGrid(transform(geometry, 3081), 5000),4269) as geometry

from tx_county

where name = "Smith";

The "uncovered" section of the county (in purple) is shown in this screenshoot with the blue arrow.  At that point, the uncovered section is right at 1000 meters north-south.

I can certainly work around the bug and get what I need, so no issue on my part now that I have a workaround.

Thanks for you help!

[sandro furieri]

I can definitely confirm: there were two different bugs affecting both
Hexagonal and Triangular grids.

A) all negative coordinate values were incorrectly handled when using
   the default origin (0,0)
   ... sorry for our spanish, north american, african and australian friends.
   latino-americans were obviously in a doubly-unlucky desperate condition:-D

B) the loops generating candidate Hexagons/Triangles to be checked for a real
   intersection with the input geom sometimes started from an incorrectly
   offset position, thus generating "blind zones"

now fixed and immediately available from the libspatialite Fossil repository

thanks for reporting this issue,
Sandro

[Shane Woodard]

Thank you Sandro!