2

I am attempting to draw a 3 km ring around downtown Minneapolis. Let downtown be defined as Government Plaza (-93.2661, 44.9765). I am using this ring to calculate census tracts within a distance of a point, so I am using the Census Bureau's projection; the census shapes are given in lat/lon coordinates, so I have to convert back and forth to calculate shape intersects on one hand, and measure distance on the other.

In Python, using shapely and pyproj:

import shapely.geometry as shpgeo, pyproj

lon, lat = (-93.2661, 44.9765)

proj = pyproj.Proj(proj="utm", ellps="GRS80", datum="NAD83")

# project into UTM to calculate distance in meters
point = shpgeo.Point(*proj(lon, lat))

# get a radius 3000 meter circle
bufferShape = point.buffer(3000)

# convert from utm back to lon/lat
coords = shp.__geo_interface__['coordinates']
# this works if shp.__geo_interface__['type'] == Polygon
newCoords = [[proj(*point, inverse=inv) for point in linring] for linring in coords] 
bufferShape = shpgeo.shape({'type': shpType, 'coordinates': tuple(newCoord)})

Now, let us test this shape, to see if we got what we wanted:

print(bufferShape.centroid)
print(bufferShape.bounds)

POINT (-93.26610002878017 44.97650391330887)
(-93.29303779998378, 44.95738308641549, -93.23916202385492, 44.99562317086678)

The centroid is just about dead on through the conversion to and from UTM; at least accurate enough for my purposes. The bufferShape should be a circle around the starting point; the max_x coord and the centroid's y coord should be the north-most point of the circle. That point is (44.9956, -93.2661); which according to distance checks with third party software is 2.12 km away; not 3 km.

How did this happen? The number 2.12 looked suspicous to me, so I changd this line:

 bufferShape = point.buffer(3000 * 1.41421)

 POINT (-93.26610005756029 44.97650782661781)
 (-93.30419574627818, 44.94946643426645, -93.22800390139928, 45.003546080298214)

The new north-most point is (-93.2661, 45.0035). Multiplying the distance by the square root of 2 in UTM, I get the correct coordinates in lon/lat.

Why is this happening?

2

The most likely answer here is that the UTM initialization is incorrect. Universal Transverse Mercator is a family of projections and by not specifying the correct zone the pyproj assumes central meridian at -183°:

Proj('+proj=tmerc +lat_0=0 +lon_0=-183 +k=0.9996 +x_0=500000 +y_0=0 +datum=NAD83
 +units=m +no_defs', preserve_units=True)

Which is likely why the grid coordinates and buffer are incorrect - this central meridian is a long way from where you are working. Instead, try specifying the correct UTM zone as part of the proj call:

proj = pyproj.Proj(proj="utm", zone=15, ellps="GRS80", datum="NAD83")

or

proj = pyproj.Proj("+init=EPSG:26915")
  • This is the correct answer. I tested in Seattle and Miami with their respective UTM zones (10 and 17) to cover the range of the continental US. The challenge now is to use the correct UTM zone for each region I am analyzing. – kingledion Jun 26 at 23:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.