Search code examples
locationgislatitude-longitudepostgisgeo

Retrieve all latitude longitude points within x distance given an area

Given an area defined by x amount of points, I'd like to split that area in to 3 mile blocks and retrieve the center latitude longitude from each block. Here's what I mean:

Area covering London:

Area split in to 3 mile blocks (not perfect):

And then print out the latitude longitude of each blocks center. This will give me all latitude longitude points (within 3mi) that cover the London area.

I need to do this programmatically, given any area of any size but unfortunately I have no idea where to start. Math isn't my strongest subject and I've not done much with geo. I think Universal Transverse Mercator coordinate system can help me, but again I don't know where to start.

Solution

  • You can actually do this as a single SQL query in Postgis, namely:

    SELECT 
  ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(
     500000 + x * 4828.03, 155000 + y * 4828.03), 27700), 4326)) 
FROM 
   generate_series(0, 12) x, 
   generate_series(0, 10) y;

    Some explanation is in order:

    1. Use Postgres's very useful generate_series function to create a grid in both the x and y directions
    2. Pass the x and y values created to ST_MakePoint to create the grid. I have used the British National Grid, which is projected and in meters, to aid the calculation. The bottom, left corner of a bounding box containing the M25 is roughly (500000, 155000), so I have used that as the anchor point. The upper right is approximately (560000, 205000), ie, 60 x 50 kilometers. Three miles is 4828.03 meters, which explains the offset in each iteration of ST_MakePoint. The limits to generate_series are 12 and 10, or 13 and 11 increments starting from 0, which is the number of times 4828.03 meters fits into the width/height of the M25's bounding box.
    3. Use ST_SetSRID to tell Postgis that these coordinates are in British National Grid, which is known as 27700.
    4. Use ST_Transform to convert to lat/lon, which is known by it's spatial reference id (SRID) as 4326.
    5. ST_AsText gives you a textual representation of the underlying geometries. Remove this to get raw geometries, ie, if you actually want to make a table and store the points as geometries.

    I used the British National Grid (27700) to create the 3 mile grid, and then converted to lat/lon. You could do the above directly by calculating roughly what 3 miles is in lat and lon in London and doing the same generate_series twice procedure, but using 4326 directly, and thereby skipping ST_Transform. The problem with this is that a minute of latitude will change as you go North/South. The British national grid has been designed with an underlying model of the earth's shape, (geoid), that is particularly well suited to the UK, and is not the same as WGS84, as used in GPS satellites. So, the logic was that using the grid in meters and then converting to lat/lon would produce more regular distortion than using lat/lon directly, where the scale changes as you move North/South. In practice, in an area the size of London, you can probably ignore all of this, unless you are building a nuclear power station :D

    You mention Spherical Mercator (3857), which is the projection used by Google Maps (and others) for global coverage in meters. You could use the above procedure using 3857 to 4326 also, but this projection maintains the increasing scale distortion as you go Northwards, whereas the UK's grid, being based on Transverse Mercator, has constant scale North/South, and distortion growing as you go East/West. If you want to plot the points on Google Maps, or similar, this might be appropriate for you. I leave it as an exercise, to figure out the start point in 3857 (hint use ST_Transform).