1

I'm trying to 1) define a square-ish boundary on a map and 2) divide that shape into a grid consisting of 1-square-mile chunks. I'm doing this because I have a dataset of people's lat/long coordinates and another dataset of business lat/long coords and am looking to simplify the calculation of distances to certain businesses for each individual (so, grouping individuals into 1-square-mile grids as opposed to treating each one individually).

Regarding 1), I've defined the four points creating my "square" below: nw_point, sw_point, ne_point, and se_point.

Regarding 2), I start at nw_point and have been trying to increment latitude and longitude appropriately to make each grid piece. I'm using a formula from this answer to build a square bounding box around each lat/long point, but am running into issues around certain latitudes. See the results corresponding to box_id 45 at the bottom – the calculated longitudinal values seem much higher than the points/boxes preceding, and this is confirmed visually, as I've been plotting the points using this website. So I'm not sure if I'm just mis-using the website, misunderstanding the type of projection being used, or whether I'm going about this grid calculation in a technically wrong way.

import numpy as np

def add_latitude(lat, mi = 1):
    modifier = mi / 69
    return lat - modifier, lat + modifier


def add_longitude(lat, long, mi = 1):
    modifier = mi / 69 / np.cos(lat)
    return long - modifier, long + modifier


def bounding_coords(lat, long, mi = 1):
    southernmost_lat, northernmost_lat = add_latitude(lat, mi = mi)
    westernmost_long, easternmost_long = add_longitude(lat, long, mi = mi)
    sw_point = (southernmost_lat, westernmost_long)
    nw_point = (northernmost_lat, westernmost_long)
    ne_point = (northernmost_lat, easternmost_long)
    se_point = (southernmost_lat, easternmost_long)
    return sw_point, nw_point, ne_point, se_point


def format_points_for_website(points, color = 'red', label = ''):
    return '\n'.join(f'{p[0]},{p[1]},{color},marker,"{label}"' for p in points)


Start with handling the contiguous United States,


by picking maximal lat/longs that completely bound the lower 48


https://en.wikipedia.org/wiki/List_of_extreme_points_of_the_United_States


Northwest Angle Inlet, MN


continental_northernmost_lat = 49.38293539482664


Ballast Key, FL


continental_southernmost_lat = 24.52108687199902


Bodelteh Islands, WA


continental_westernmost_long = -124.76410018717496


Sail Rock, Lubec, ME


continental_easternmost_long = -66.94700844863216
nw_point = (continental_northernmost_lat, continental_westernmost_long)
sw_point = (continental_southernmost_lat, continental_westernmost_long)
ne_point = (continental_northernmost_lat, continental_easternmost_long)
se_point = (continental_southernmost_lat, continental_easternmost_long)


Starting at northwesternmost point above, begin building out "chunks" of land.


You could technically start at any of the four points, but we'll start with this one.


Start by traversing land until we've passed the southernmost point, then increment to the east


and start all over again until we've traversed to the easternmost point.


curr_long = nw_point[1]
begin_lat, begin_long = nw_point
boxes = {}
box_id = 0
colors = ['red', 'green', 'blue', 'purple', 'orange']
points = []


Until we've passed the easternmost continental longitude


while curr_long < continental_easternmost_long:
    # And until we've passed the southernmost latitude
    curr_lat = begin_lat
    while curr_lat > continental_southernmost_lat:
        if box_id == 46:raise
        # First or 0th element represents a numerically lower latitude, which we need to
        # use, since latitude increases the further north you go, and we're building from northwest to southeast
        curr_lat = add_latitude(curr_lat)[0]
        # Get lower latitude boundary for use in centroid of "next" bounding box
        lower_lat = add_latitude(curr_lat)[0]
        print(f'Center point: ({lower_lat}, {curr_long})')
        sw_point, nw_point, ne_point, se_point = bounding_coords(lower_lat, curr_long)
        print(format_points_for_website([sw_point, nw_point, ne_point, se_point], color = colors[box_id % len(colors)], label = box_id))
        # Could use a more informative box ID to perhaps represent common lat/long
        # boundaries across different boxes, but this suffices for now...
        box_id += 1
    # Use second element (numerically greater)
    curr_long = add_longitude(curr_lat, curr_long)[1]

Relevant output provided below:

...
Center point: (48.7307614817832, -124.76410018717496)
48.71626872816001,-125.16592296075729,purple,marker,"43"
48.745254235406385,-125.16592296075729,purple,marker,"43"
48.745254235406385,-124.36227741359264,purple,marker,"43"
48.71626872816001,-124.36227741359264,purple,marker,"43"
Center point: (48.71626872816001, -124.76410018717496)
48.70177597453682,-125.43565414695419,orange,marker,"44"
48.7307614817832,-125.43565414695419,orange,marker,"44"
48.7307614817832,-124.09254622739573,orange,marker,"44"
48.70177597453682,-124.09254622739573,orange,marker,"44"
Center point: (48.70177597453682, -124.76410018717496)
48.687283220913635,-126.80827425555613,red,marker,"45"   ***
48.71626872816001,-126.80827425555613,red,marker,"45"   ***
48.71626872816001,-122.7199261187938,red,marker,"45"   ***
48.687283220913635,-122.7199261187938,red,marker,"45"   ***
boot-scootin
  • 111
  • 4
  • 2
    Calculating true distance is going to be easier than trying to compute a "square" where the northern pair of points are closer together than the southern pair (in the Northern Hemisphere; reversed in the Southern Hemisphere). While it isn't impossible to generate a rectangular grid of nearly equal area sections, you have to choose one dimension to be fixed interval, and that pretty much has to be longitude, and you have to be open to the fact the objects are spheroidal trapezoids, not rectangles. This seems to be an XY Problem. – Vince Sep 02 '21 at 14:25
  • Would holding longitude on a fixed interval and varying latitude guarantee that all area in my "square" is assigned to one grid section? Or could there be gaps where some areas are left out? – boot-scootin Sep 02 '21 at 14:39
  • If you partition a space in 2D, then it's partitioned by 2D functions. It may not be partitioned by 3D functions. Note that plotting {lat,lon} values is plotting {Y,X}, which may give some mapping apps fits. – Vince Sep 02 '21 at 14:46
  • 1
    Maybe it would help to use shapely geometries in combination with geopandas, so you can project points in a certain region to a metric CRS. Once you have that, you can easily obtain square bounding boxes and divide them into equal parts. – Ben Sep 03 '21 at 10:41
  • @Ben why would projecting my points into a metric CRS solve the problem, vs. a degree-based CRS? – boot-scootin Sep 03 '21 at 13:22

1 Answers1

1

I am going to put this as an answer because I'll post some code.

@blacksite (in the comments): It may be easier to work in a metric crs. Consider this, where I convert your point to UTM zone 10N (EPSG:32610), to get a metric projection, then create a circular buffer and take its bounding box:

import geopandas as gpd
from shapely.geometry import Point, box

p = Point((-124.76410018717496, 48.7307614817832))  # lon, lat
pp = gpd.GeoSeries(p).set_crs(4326).to_crs(32610)[0] # Convert to UTM for plot)


gs = gpd.GeoSeries(p).set_crs(4326).to_crs(32610).buffer(1609,34) # 1 mile circular buffer


gs_plot = gpd.GeoSeries([pp,box(*gs.bounds.values[0])]).set_crs(32610) # For plotting


Plot


ax = gs_plot[1:].plot()
gs_plot[:1].plot(ax=ax, color="red")

enter image description here

ax = gs_plot.to_crs(4326)[1:].plot()
gs_plot[:1].to_crs(4326).plot(ax=ax, color="red")

enter image description here

I think, this is much easier. The only "difficult" part now would be to check for each of your points, in which UTM zone they lie in, but I'm pretty sure that there's a very nice way to find out by the lat/lon coordinates.

I hope this helps. I am not entirely sure whether this is what you are looking for.

Ben
  • 193
  • 1
  • 8