Stochastic Coder

I was working on a Python project in healthcare. I needed to find gaps in emergency transportation coverage based on specific locations and service ranges. The coverage could vary for each healthcare provider and mode of transport (air rotor, air fixed wing, land vehicle) combination. I wanted to represent the coverage with circles, where the radius indicates coverage in miles. However, I also needed control over certain rendering elements. To complete the project, I had to create and display coverage circles for each location, so I needed to create a function for this. The function for rendering each circle is based on the Haversine Formula.

Example using Haversine Formula in Python

For this example, I choose Memorial Stadium, the venue of the Nebraska Huskers, my favorite college football team, situated in the heart of the United States.

The map above shows a Basemap using an ‘orthographic’ projection (ortho for short) with a 50-mile radius drawn around Memorial Stadium. This helps you to understand the location of Memorial Stadium and Nebraska (the state within the green radius).

Let’s look at the code

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import math
from osgeo import ogr, osr
import matplotlib.cm as cmx
import matplotlib.colors as colors

The Haversine Functions

def createCircleAroundWithRadius(lat, lon, radiusMiles):
 ring = ogr.Geometry(ogr.wkbLinearRing)
 latArray = []
 lonArray = []
 
 for brng in range(0,360):
   lat2, lon2 = getLocation(lat,lon,brng,radiusMiles)
   latArray.append(lat2)
   lonArray.append(lon2)

return lonArray,latArray


def getLocation(lat1, lon1, brng, distanceMiles):
 lat1 = lat1 * math.pi/ 180.0
 lon1 = lon1 * math.pi / 180.0
 #earth radius
 #R = 6378.1Km
 #R = ~ 3959 MilesR = 3959

 distanceMiles = distanceMiles/R

 brng = (brng / 90)* math.pi / 2

 lat2 = math.asin(math.sin(lat1) * math.cos(distanceMiles) 
   + math.cos(lat1) * math.sin(distanceMiles) * math.cos(brng))

 lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(distanceMiles)
   * math.cos(lat1),math.cos(distanceMiles)-math.sin(lat1)*math.sin(lat2))

 lon2 = 180.0 * lon2/ math.pi
 lat2 = 180.0 * lat2/ math.pi

return lat2, lon2

The function createCircleAroundWithRadius takes the latitude and longitude of the circle’s center, as well as the distance from the center at which the radius should be drawn (which can also be specified in kilometers). This function primarily iterates over bearings from 0 to 360 degrees and determines the radius points within the getLocation function.

Orthographic Projection

memorialStaduimLat = 40.820485
memorialStadiumLon = -96.705588
distanceInMiles = 50

#orthographic basemap, center on specified lat and lon
mp = Basemap(projection='ortho', lat_0 = 44, lon_0 = -105, resolution = 'l')


mp.drawstates()
mp.drawcountries()
mp.drawmapboundary(fill_color='aqua')
mp.fillcontinents(color='coral',lake_color='aqua')
mp.drawcoastlines()

#color radius green - seems like a fine color
c = 'g'
# retrieve X and Y radius values using memorial stadium as center point, draw 50 miles out
X,Y = createCircleAroundWithRadius(memorialStaduimLat,memorialStadiumLon,
distanceInMiles)

X,Y = mp(X,Y)
mp.plot(X,Y,marker=None,color=c,linewidth=2)

x,y = mp(memorialStadiumLon,memorialStaduimLat)
mp.plot(x,y ,marker='D',color=c,markersize=4)

Mercator Projection

Most of the maps I create don’t use the orthographic projection, but instead use the Mercator (merc) projection. Starting with the map below might have caused confusion among some readers about the specific location shown. Orthographic projections can be helpful when points are spread across a wide area, or when a particular point is unfamiliar to the audience (as in this case). The geospatial reports I create usually focus on states and counties, which is why I often choose the Mercator projection. Others may have had different experiences.

westbounds = -105.053514
eastbounds = -95.30829
northbounds = 44.001707
southbounds = 39.999932fig = plt.figure(1)

#merc basemap, center on specified lat and lon, and zoom to bounds
mp = Basemap(projection='merc', lat_0 = 44, lon_0 = -105, ax=fig.gca(),
resolution = 'h', area_thresh = 0.1,
llcrnrlon=westbounds, llcrnrlat=southbounds,
urcrnrlon=eastbounds, urcrnrlat=northbounds)

mp.drawstates()
mp.drawcounties()
memorialStaduimLat = 40.820485
memorialStadiumLon = -96.705588
distanceInMiles = 50

#color radius green - seems like a fine color
c = 'g'

# retrieve X and Y radius values using memorial stadium as center point, draw 50 miles out
X,Y = createCircleAroundWithRadius(memorialStaduimLat,memorialStadiumLon,distanceInMiles)
X,Y = mp(X,Y)
mp.plot(X,Y,marker=None,color=c,linewidth=2)

x,y = mp(memorialStadiumLon,memorialStaduimLat)
mp.plot(x,y ,marker='D',color=c,markersize=4)

I trust this will assist anyone seeking to denote coverage areas or distances using circles and Basemap. Please do not hesitate to reach out if you have any inquiries or issues.