How can I generate a regular geographic grid using python?

Question

I want to retrieve all lat/lon coordinate pairs of a regular grid over a certain map area. I have found the geopy library, but didn't manage at all to approach the problem.

For example, I have a rectangular geographic area described by its four corners in lat/lon coordinates, I seek to calculate the grid with spacing of e.g. 1km covering this area.

Spit

Cardinal points like NSEW? How can cardinal points be given in (lon, lat)-pairs and how can they have a distance? Please elaborate more on your problem and add the code you have so far. – Hippocras

Sorry, i was wrong i want to find all lat long coordinates in a certain area, for example i have 4 points lat long of a region in google maps, i want to retrive all coordinates lat long in this area with 10 kilometers of distance each other. i dont have a code now, im still searching something in geopy. – Spit

So you want to create a regular grid of point coordinates given a certain polygon? – Hippocras

yes, is this! is there anything in geopy or in other library ? – Spit

Answer 1

Preliminary Considerations

It makes a slight difference how your certain area is defined. If it's just a rectangular area (note: rectangular in projection is not neccessarily rectangular on Earth's surface!), you can simply iterate from min to max in both coordinate dimensions using your desired step size. If you have an arbitrary polygon shape at hand, you need to test which of your generated points intersects with that poly and only return those coordinate pairs for which this condition holds true.

Calculating a Regular Grid

A regular grid does not equal a regular grid across projections. You are talking about latitude/longitude pairs, which is a polar coordinate system measured in degrees on an approximation of Earth's surface shape. In lat/lon (EPSG:4326), distances are not measured in meters/kilometers/miles, but rather in degrees.

Furthermore, I assume you want to calculate a grid with its "horizontal" steps being parallel to the equator (i.e. latitudes). For other grids (for example rotated rectangular grids, verticals being parallel to longitudes, etc.) you need to spend more effort transforming your shapes.

Ask yourself: Do you want to create a regularly spaced grid in degrees or in meters?

A grid in degrees

If you want it in degrees, you can simply iterate:

stepsize = 0.001
for x in range(lonmin, lonmax, stepsize):
    for y in range(latmin, latmax, stepsize):
        yield (x, y)

But: Be sure to know that the length in meters of a step in degrees is not the same across Earth's surface. For example, 0.001 delta degrees of latitude close to the equator covers a different distance in meters on the surface than it does close to the poles.

A grid in meters

If you want it in meters, you need to project your lat/lon boundaries of your input area (your certain area on a map) into a coordinate system that supports distances in meters. You can use the Haversine formula as a rough approximation to calculate the distance between lat/lon pairs, but this is not the best method you can use.

Better yet is to search for a suitable projection, transform your area of interest into that projection, create a grid by straightforward iteration, get the points, and project them back to lat/lon pairs. For example, a suitable projection for Europe would be EPSG:3035. Google Maps uses EPSG:900913 for their web mapping service, by the way.

In python, you can use the libraries shapely and pyproj to work on geographic shapes and projections:

import shapely.geometry
import pyproj

# Set up transformers, EPSG:3857 is metric, same as EPSG:900913
to_proxy_transformer = pyproj.Transformer.from_crs('epsg:4326', 'epsg:3857')
to_original_transformer = pyproj.Transformer.from_crs('epsg:3857', 'epsg:4326')

# Create corners of rectangle to be transformed to a grid
sw = shapely.geometry.Point((-5.0, 40.0))
ne = shapely.geometry.Point((-4.0, 41.0))

stepsize = 5000 # 5 km grid step size

# Project corners to target projection
transformed_sw = to_proxy_transformer.transform(sw.x, sw.y) # Transform NW point to 3857
transformed_ne = to_proxy_transformer.transform(ne.x, ne.y) # .. same for SE

# Iterate over 2D area
gridpoints = []
x = transformed_sw[0]
while x < transformed_ne[0]:
    y = transformed_sw[1]
    while y < transformed_ne[1]:
        p = shapely.geometry.Point(to_original_transformer.transform(x, y))
        gridpoints.append(p)
        y += stepsize
    x += stepsize

with open('testout.csv', 'wb') as of:
    of.write('lon;lat\n')
    for p in gridpoints:
        of.write('{:f};{:f}\n'.format(p.x, p.y))

This example generates this evenly-spaced grid:

Hippocras

thanks a lot, you have been very clear, and you taught me many new things, thank you! – Spit

There is a little error in your code: when writing out to testout.csv, the stream is defined as 'wb', but a string is written. It should be changed to 'w'. – Parlance

Right, it can be changed to w here; however, since this is Python 2.7 I, personally, usually go with wb when writing strings to correctly write Unicode. – Hippocras

Great stuff jbndlr - I am trying to build something where a user can define this area themselves and therefore I need a solution to use the correct projections in an automated way - any ideas? – Saunders

Depends on your actual needs; did you take a look at the QGIS toolboxes? If your users draw their polygons on a map, then the projection should be clear; if not, you need to handle different srids depending on input projection. QGIS (et al.) can surely help you! – Hippocras

How would you go about creating patches to then check if a certain coordinate lands in a patch? – Stutz

Same approach if you want to do it manually (make sure grid and your coordinate have the same projection) - just instead of creating Point shapes, create Polygon shapes with +- 1/2 distance in all four directions. Or use something like Google's S2 or Uber's H3. – Hippocras

you mixed up the corner names. nw is not the top left corner but rather the bottom left corner and therefore sw. Same for se which should be named ne for north-east corner (top right). – Same

How are you plotting this visualization? The only output I am seeing is the csv you are creating? – Tyrannize

@Tyrannize I loaded the resulting CSV in QGIS and took a screenshot, that was the quickest solution for me. You can also easily use matplotlib et al – Hippocras

This seems to not work for really small areas where the stepsize is small. I tried this with a much smaller sw and ne point with a smaller step size but it produced zero results. Any thoughts? – Tyrannize

Hum, maybe precision limitations on Lon/Lat in combination with PROJ.4's reprojection? – Hippocras

Great solution. There is a minor mistake though - to_original_transformer = pyproj.Transformer.from_crs('epsg:4326', 'epsg:3857') Should be: to_original_transformer = pyproj.Transformer.from_crs('epsg:3857', 'epsg:4326') – Fidole

I have a question. The code as presented above works fine. But with bounds of New York: sw = shapely.geometry.Point((-74.25, 40.49)) ne = shapely.geometry.Point((-73.65, 40.92)) It does not work properly and gives me distances of 1.35km instead of 5km. Why is that? – Buckles