I was working on a project were I needed to identify gaps in coverage based on a location and a service range, which could be different for each service/provider combination. A coverage circle where the radius represents coverage miles would suffice for this project, but I required control over some rendering elements. To complete the project I would have to draw my own coverage circles for each location, thus I needed a function. The function to render each circle is based on the Haversine Formula (https://en.wikipedia.org/wiki/Haversine_formula) if anyone is ever really interested in the steps please let me know and I will include additional details. Otherwise, I will just stick with the code.
50 Miles From Memorial Stadium, Lincoln NE, USA
Well I could not think of a really good well known ‘neutral’ location to use in my example, so I picked the most singular location where I live, Memorial Stadium, home of the UNL football team; ‘Smack dab’ in the middle of the U.S.
Above is a Basemap plot using an ‘orthographic’ project (ortho for short) with a 50 mile radius drawn around Memorial Stadium. If nothing else, you now know where Memorial Stadium is and Nebraska is (the state containing 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 expects the lat and long of the center of the circle as well as the miles from the center the radius should be drawn (this could easily be kilometers). This function does little else except iterate over barrings 0-360 degrees. The determination of the radius points is actually done in the function getLocation.
Let’s draw something – 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)
How about a Mercator Projection
Most of the maps I render do not use the ortho project, but instead use the Mercator (merc) projection. If I would have started with the map below some readers may have been confused as to where in the world we are looking. Orthographic projections can be great when points are spread across a very large area, or when a relevant point is not familiar to an audience (as in this case). Usually, geospatial reports I produce focus on states and counties which is why I often use 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 hope this help anyone wanting to indicate coverage areas or distances using circles and Basemap. I am thinking about doing something similar in R if need be. Please let me know if you have any questions or concerns.