LonLat Points to Grid of Polygons in PostGIS

Question

I'm going to preface this with it's my first time using PostGIS.

I'm starting out with a 5 km x 5 km grid cell in which I know the center point lon/lat (SRID=4326). I'm trying to create a polygon for each of these grid cells in PostGIS using this command:

    create table lonlat_poly as
    select ST_MakeEnvelope(lon-2.5 km, lat-2.5 km, lon+2.5 km, lat+2.5 km)
    as newlonlat from lonlat;

where lonlat contains the two columns lon and lat.

The potential pitfall is that the simple subtraction and addition of 2.5 km is not going to be correct due to projection issues (right?). I'm sure there's a more precise way of doing this, but I haven't thought or stumbled upon anything yet. Thanks in advance for any hints you can provide!

score 5 · Answer 1 · edited Apr 13 '17 at 12:34

You have two options: make a perfect 5 km grid, projected inperfectly to lat/long, or make a perfect lat/long grid that approximates a 5 km grid. See How to create a regular polygon grid in PostGIS? to help make a grid.

Perfect 5 km grid

For example, using this function, here is a 4 x 6 grid with 5 km x 5 km cells somewhere in Vancouver using NAD83 UTM Zone 10N (EPSG:26910), then transformed back into WGS 84 lat/long (EPSG:4326):

SELECT ST_Transform(ST_SetSRID(ST_Collect(cells), 26910), 4326)
FROM ST_CreateFishnet(4, 6, 5000, 5000, 493000, 5456000) AS cells;

You can approximate your x0 and y0 origin reference from your projected grid points (subtracting 2500 m from each x and y).

Perfect lat/long grid

Similar to above, but you must approximate 5 km in degrees latitude and longitude. (I just used one point, but you might want to approximate this distance at the centroid of your points).

SELECT ST_SetSRID(ST_Collect(cells), 4326)
FROM ST_CreateFishnet(4, 6, 0.0687135352095627, 0.0449759180267531, -123.096198961667, 49.2567270585716) AS cells;

Both results are pretty close to each other due to the small scale, but will look more warped over larger regions.

Thanks @Mike Toews I like these approaches. My Grid covers the whole US, so I opted to use the ST_Transform function on each lon lat pair and then add and subtract my 2.5 km from that transform and it looks pretty good. — Jonathan, May 10 '12 at 17:41

score 2 · Accepted Answer · answered May 09 '12 at 19:31

2

you are right, sometimes projections cant make result as we want. in this case, i âm converting my projection to Google mercator for getting true result.

i dont know is it work for you but i use ST_Transform method like this case.

    create table lonlat_poly as
    select ST_Transform(ST_MakeEnvelope(lon-2.5 km, lat-2.5 km, lon+2.5 km, 
    lat+2.5 km),900913)
    as newlonlat from lonlat;

i hope it helps you....

answered May 09 '12 at 19:31

urcm

22,533
4
57
109

Thanks @Aragon: I used this approach, but also had to add an additional ST_Transform within the MakeEnvelope command to convert each of the lon and lat values to meters. Thanks again! – Jonathan May 10 '12 at 17:38

LonLat Points to Grid of Polygons in PostGIS

2 Answers2