Designing a grid overlay based on longitudes and latitudes

I'm trying to figure out the best way to approach the following:

Say I have a flat representation of the earth. I would like to create a grid that overlays this with each square on the grid corresponding to about 3 square kilometers. Each square would have a unique region id. This grid would just be stored in a database table that would have a region id and then probably the long/lat coordinates of the four corners of the region, right? Any suggestions on how to generate this table easily? I know I would first need to find out the width and height of this "flattened earth" in kms, calculate the number of regions, and then somehow assign the long/lats to each intersection of vertical/horizontal line; however, this sounds like a lot of manual work.

Secondly, o开发者_StackOverflow中文版nce I have that grid table created, I need to design a fxn that takes a long/lat pair and then determines which logical "region" it is in. I'm not sure how to go about this.

Any help would be appreciated.

Thanks.

---

Assume the Earth is a sphere with radius R = 6371 km.

Start at (lat, long) = (0, 0) deg. Around the equator, 3km corresponds to a change in longitude of

dlong = 3 / (2 * pi * R) * 360 
      = 0.0269796482 degrees

If we walk around the equator and put a marker every 3km, there will be about (2 * pi * R) / 3 = 13343.3912 of them. "About" because it's your decision how to handle the extra 0.3912.

From (0, 0), we walk North 3 km to (lat, long) (0.0269796482, 0). We will walk around the Earth again on a path that is locally parallel to the first path we walked. Because it is a little closer to the N Pole, the radius of this circle is a bit smaller than that of the first circle we walked. Let's use lower case r for this radius

r = R * cos(lat)
  = 6371 * cos(0.0269796482)
  = 6 368.68141 km

We calculate dlong again using the smaller radius,

dlong = 3 / (2 * pi * r) * 360
      = 0.0269894704 deg

We put down the second set of flags. This time there are about (2 * pi * r) / 3 = 13 338.5352 of them. There were 13,343 before, but now there are 13,338. What's that? five less.

How do we draw a ribbon of squares when there are five less corners in the top line? In fact, as we walked around the Earth, we'd find that we started off with pretty good squares, but that the shape of the regions sheared out into pretty extreme parallelograms.

We need a different strategy that gives us the same number of corners above and below. If the lower boundary (SW-SE) is 3 km long, then the top should be a little shorter, to make a ribbon of trapeziums.

There are many ways to craft a compromise that approximates your ideal square grid. This wikipedia article on map projections that preserve a metric property, links to several dozen such strategies.

The specifics of your app may allow you to simplify things considerably, especially if you don't really need to map the entire globe.

---

Microsoft has been investing in spatial data types in their SQL Server 2008 offering. It could help you out here. Because it has data types to represent your flattened earth regions, operators to determine when a set of coordinates is inside a geometry, etc. Even if you choose not to use this, consider checking out the following links. The second one in particular has a lot of good background information on the problem and a discussion on some of the industry standard data formats for spatial data.

http://www.microsoft.com/sqlserver/2008/en/us/spatial-data.aspx

http://jasonfollas.com/blog/archive/2008/03/14/sql-server-2008-spatial-data-part-1.aspx

---

First, Paul is right. Unfortunately the earth is round which really complicates the heck out of this stuff.

I created a grid similar to this for a topographical mapping server many years ago. I just recoreded the coordinates of the upper left coder of each region. I also used UTM coordinates instead of lat/long. If you know that each region covers 3 square kilometers and since UTM is based on meters, it is straight forward to do a range query to discover the right region.

---

You do realize that because the earth is a sphere that "3 square km" is going to be a different number of degrees near the poles than near the equator, right? And that at the top and bottom of the map your grid squares will actually represent pie-shaped parts of the world, right?

I've done something similar with my database - I've broken it up into quad cells. So what I did was divide the earth into four quarters (-180,-90)-(0,0), (-180,0)-(0,90) and so on. As I added point entities to my database, if the "cell" got more than X entries, I split the cell into 4. That means that in areas of the world with lots of point entities, I have a lot of quad cells, but in other parts of the world I have very few.

My database for the quad tree looks like:

\d areaids;
                 Table "public.areaids"
    Column    |            Type             | Modifiers 
--------------+-----------------------------+-----------
 areaid       | integer                     | not null
 supercededon | timestamp without time zone | 
 supercedes   | integer                     | 
 numpoints    | integer                     | not null
 rectangle    | geometry                    | 
Indexes:
    "areaids_pk" PRIMARY KEY, btree (areaid)
    "areaids_rect_idx" gist (rectangle)
Check constraints:
    "enforce_dims_rectangle" CHECK (ndims(rectangle) = 2)
    "enforce_geotype_rectangle" CHECK (geometrytype(rectangle) = 'POLYGON'::text OR rectangle IS NULL)
    "enforce_srid_rectangle" CHECK (srid(rectangle) = 4326)

I'm using PostGIS to help find points in a cell. If I look at a cell, I can tell if it's been split because supercededon is not null. I can find its children by looking for ones that have supercedes equal to its id. And I can dig down from top to bottom until I find the ones that cover the area I'm concerned about by looking for ones with supercedeson null and whose rectangle overlaps my area of interest (using the PostGIS '&' operator).

---

There's no way you'll be able to do this with rectangular cells, but I've just finished an R package dggridR which would make this easy to do using a grid of hexagonal cells. However, the 3km cell requirement might yield so many cells as to overload your machine.

You can use R to generate the grid:

install.packages('devtools')
install.packages('rgdal')
library(devtools)
devools.install_github('r-barnes/dggridR')
library(dggridR)
library(rgdal)

#Construct a discrete global grid (geodesic) with cells of ~3 km^2
dggs <- dgconstruct(area=100000, metric=FALSE, resround='nearest')

#Get a hexagonal grid for the whole earth based on this dggs
grid <- dgearthgrid(dggs,frame=FALSE)

#Save the grid
writeOGR(grid, "grid_3km_cells.kml", "cells", "KML")

The KML file then contains the ids and edge vertex coordinates of every cell.

The grid looks a little like this:

My package is based on Kevin Sahr's DGGRID which can generate this same grid to KML directly, though you'll need to figure out how to compile it yourself.