2

I was reviewing this function in PostGis

https://postgis.net/docs/manual-dev/ST_HexagonGrid.html

1) What I don't understand is what the underlying geom data would be. What's the source to get the USA map as shown? What is the DB schema? I think it could be one record, if I only need the USA boundary, not each state?

2) Is the result a list of points? or geom vectors?

3) If geom vectors, how do you convert them into points of lat and lng?

4) How to do approximate the hexigons to approximate a 50 mile radius from a point?

UPDATE:

I played with the width to try to get the correct number of hexagons based on Jim Jones example below. Unfortunately, something went wrong..

1) the length seems to have no relation to meters

2) There are multiple sized hexagons, which seems weird.

postgis_test=# WITH j AS (
postgis_test(# SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
postgis_test(#   SELECT 
postgis_test(#   generate_hexgrid(
postgis_test(#     5909968.8,
postgis_test(#     ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
postgis_test(# FROM usa_states)i) 
postgis_test-# SELECT count(j.hex) FROM j,usa_states
postgis_test-# WHERE ST_Intersects(usa_states.geom,j.hex);
 count 
-------
   119
(1 row)

postgis_test=# WITH j AS (
postgis_test(# SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
postgis_test(#   SELECT 
postgis_test(#   generate_hexgrid(
postgis_test(#     5909968.8,
postgis_test(#     ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
postgis_test(# FROM usa_states)i) 
postgis_test-# SELECT DISTINCT st_area(j.hex) FROM j,usa_states
postgis_test-# WHERE ST_Intersects(usa_states.geom,j.hex);
     st_area      
------------------
 1219.78281686003
 2089.11341619338
 2089.11341619338
 3379.93344444246
  7051.4650344734
 12076.9943663072
(6 rows)
user2012677
  • 5,465
  • 6
  • 51
  • 113
  • Some time ago I adapted a function to create hexagons to answer this question: https://stackoverflow.com/a/50845811/2275388 does it help you? – Jim Jones Jan 21 '20 at 20:11
  • Great project! How would one make sure the sizes of the hexagons are equal distance/area and that they can be used a proxy for 50-mile radius? – user2012677 Jan 21 '20 at 20:20
  • Perhaps you could elaborate a bit about your use case and add to your question. Why do you need hexagons? You need hexagons that cover a certain area or you want to create an hexagon around a point - like a buffer? – Jim Jones Jan 21 '20 at 20:24
  • I have an api, and I can search by radius (50 miles). I want to cover/query the entire USA making the least number of calls to the API as possible, so I am trying to calculate what locations that I need to query. – user2012677 Jan 21 '20 at 20:28
  • so you want to cover the whole country with hexagons of the same size - 50 miles? – Jim Jones Jan 21 '20 at 20:30
  • Basically, YES! The hexigons would fit into a 50 miles radius circle, so there would be some overlap of the circles. Here is an example https://stackoverflow.com/questions/1404944/covering-an-arbitrary-area-with-circles-of-equal-radius?noredirect=1&lq=1 – user2012677 Jan 21 '20 at 20:33
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/206382/discussion-between-jim-jones-and-user2012677). – Jim Jones Jan 21 '20 at 20:45
  • I just added a function that might do what you want. I'll work on a better explanation later today. – Jim Jones Jan 22 '20 at 06:47

1 Answers1

2

According to the author, the following function should create a grid with the extent based on the given BBOX and the cell size in meters.

SRID 3857 units are [very approximately] meters, and using this projection will create hex cells that "look right" on a web map (most of which use a web mercator projection).

CREATE OR REPLACE FUNCTION generate_hexgrid(width float, xmin float, ymin float, xmax float, ymax float, srid int default 3857)
RETURNS TABLE(gid text,geom geometry(Polygon)) AS $$
DECLARE
  b float := width / 2;
  a float := tan(radians(30)) * b;
  c float := 2 * a;
  height float := 2 * (a + c);
  index_xmin int := floor(xmin / width);
  index_ymin int := floor(ymin / height);
  index_xmax int := ceil(xmax / width);
  index_ymax int := ceil(ymax / height);
  snap_xmin float := index_xmin * width;
  snap_ymin float := index_ymin * height;
  snap_xmax float := index_xmax * width;
  snap_ymax float := index_ymax * height;
  ncol int := abs(index_xmax - index_xmin);
  nrow int := abs(index_ymax - index_ymin);
  polygon_string varchar := 
    'POLYGON((' || 0 || ' ' || 0 || ' , ' || b || ' ' || a || ' , ' ||
    b || ' ' || a + c || ' , ' || 0 || ' ' || a + c + a || ' , ' ||
    -1 * b || ' ' || a + c || ' , ' || -1 * b || ' ' || a || ' , ' ||
    0 || ' ' || 0 ||'))';
BEGIN
  RETURN QUERY
  SELECT 
    format('%s %s %s', width,
    x_offset + (1 * x_series + index_xmin),
    y_offset + (2 * y_series + index_ymin)),
    ST_SetSRID(ST_Translate(two_hex.geom,
    x_series * width + snap_xmin,
    y_series * height + snap_ymin), srid)
  FROM  generate_series(0, ncol, 1) AS x_series,
        generate_series(0, nrow, 1) AS y_series,
    (SELECT 0 AS x_offset, 0 AS y_offset, polygon_string::geometry AS geom
     UNION
     SELECT 0 AS x_offset, 1 AS y_offset, ST_Translate(polygon_string::geometry, b , a + c)  AS geom
    ) AS two_hex;
END; $$ LANGUAGE plpgsql;

Considering you have a table called usa and it contains the geometries of this shapefile you should be able to create a grid that overlaps the USA map with the following query:

CREATE TABLE usa_hex AS
WITH j AS (
SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
  SELECT 
  generate_hexgrid(
    80467,
    ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
    ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
    ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
    ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
FROM usa)i) 
SELECT j.hex FROM j,usa
WHERE ST_Intersects(usa.geom,j.hex);

enter image description here

EDIT: It is still not an answer, as it does not create the hexagons using meters, but it might help other users. The following function (derived from this answer) creates geometry type hexagons of exact the same size in degrees. 

CREATE OR REPLACE FUNCTION generate_hexagons(width FLOAT, bbox BOX2D, srid INTEGER DEFAULT 4326)
RETURNS TABLE (gid INTEGER, hexagon GEOMETRY) AS $$
DECLARE
  b FLOAT := width/2;
  a FLOAT := b/2;
  c FLOAT := 2*a;
  height FLOAT := 2*a+c;
  ncol FLOAT := ceil(abs(ST_Xmax(bbox)-ST_Xmin(bbox))/width);
  nrow FLOAT := ceil(abs(ST_Ymax(bbox)-ST_Ymin(bbox))/height);
  polygon_string VARCHAR := 'POLYGON((' || 
    0 || ' ' || 0 || ' , ' || b || ' ' || a || ' , ' || b || ' ' || a+c || ' , ' || 0 || ' ' || a+c+a || ' , ' ||
   -1*b || ' ' || a+c || ' , ' || -1*b || ' ' || a || ' , ' || 0 || ' ' || 0 || '))';
BEGIN    
  RETURN QUERY 
  SELECT 
    row_number() OVER ()::INTEGER,
    ST_SetSRID(
      ST_Translate(geom, x_series*(2*a+c)+ST_Xmin(bbox), y_series*(2*(c+a))+ST_Ymin(bbox)),srid)
  FROM generate_series(0, ncol::INTEGER, 1) AS x_series,
       generate_series(0, nrow::INTEGER,1 ) AS y_series,
       (SELECT polygon_string::GEOMETRY AS geom
        UNION
        SELECT ST_Translate(polygon_string::GEOMETRY, b, a + c) AS geom) AS two_hex;    
END;
$$ LANGUAGE plpgsql;

Overlapping with the data set used above:

WITH j (hex_rec) AS (
  SELECT generate_hexagons(3.0,ST_Extent(geom)) 
  FROM usa
)
SELECT (hex_rec).gid,(hex_rec).hexagon FROM j, usa 
WHERE ST_Intersects(usa.geom,(hex_rec).hexagon);

enter image description here

Further reading:

Jim Jones
  • 18,404
  • 3
  • 35
  • 44
  • I see you hard coded a points for min x and y, and max x and y. What are those points? Also what is 3857 represent vs WGS84? – user2012677 Jan 22 '20 at 15:26
  • 1
    4326 is the number for the reference system WGS84. 3857 is another reference system :) x and y are not hard coded, they represent the extent of the BBOX. You can obviously change it to cover a smaller or bigger area :) good luck! – Jim Jones Jan 22 '20 at 18:26
  • Thank you, but I think something is wrong, there are too many hexigons. The script created over 5075 hexigons, but I estimated there should only be about 473 hexigons using a 50 mile radius, 118 using 100 miles radius. I beleive this script is using a 50mile radius. Am I miscalculating? – user2012677 Jan 23 '20 at 01:39
  • I'm unable to calculate how many 50 miles hexagons you should get from US territory right now, but have you played with the size parameter to see if it gives you a satisfactory number of hexagons? – Jim Jones Jan 23 '20 at 09:30
  • Are you recommending I just guess he width? – user2012677 Jan 23 '20 at 12:59
  • @user2012677 no, I meant you could check if the created hexagons are right and your calculations are wrong. If the hexagons are indeed too small, just increase their size by passing a bigger parameter. Take into account the SRS 3857 if you intend to change the hexagons size. Good luck. – Jim Jones Jan 23 '20 at 13:05
  • updated notes above. Seems I'm getting mutiple sized hexagons – user2012677 Jan 26 '20 at 13:02
  • @user2012677 I see the function I found wasn't that accurate after all :) If you already have found a solution by now, please post it here :) – Jim Jones Jan 30 '20 at 10:31
  • Never found solution – user2012677 Jan 30 '20 at 12:01
  • @user2012677 too bad. I'm also not aware of a function to create hexagons based on ml/km. If you can live with degrees, the first function we used (added today in the answer) does a pretty good job in creating hexagons of the same size – Jim Jones Jan 30 '20 at 12:05
  • What does st_area output to? Square meters? As I can try adjusting the hexagon in degrees, and reviewing the st_area if I know the output to compare to the expected square meters – user2012677 Jan 30 '20 at 14:22
  • @user2012677 that's the tricky part. ST_Area with `geometry` returns the area in the unit of measurement of your SRID, which in you case is `WGS84`. If your use case is only about creating a buffer to speed up your queries, it wouldn't matter that much if you created the hexagon in miles or degrees. Keep in min: If you switch to `geography` the values will be totally different. – Jim Jones Jan 30 '20 at 14:33