How to plot counties and shade them according with a specified color?
I forgot that this has been an old question for me.
Where to get county shape files?
A UScounty shapefile is included in basemap version 1.0.6. Unfortunately, the latest version of basemap provided by enpkg is 1.0.2 (and I've not successfully been able to compile from source on a Mac).
I think it's safe to copy the USCounties files over from version 1.0.6 to your basemap data directory, but I'd rather not add files to that directory.
Alternative: download files from US Census
http://www.census.gov/geo/www/tiger/tgrshp2012/tgrshp2012.html
ftp interface: ftp://ftp2.census.gov/geo/tiger/TIGER2012/COUNTY/ -> ftp://ftp2.census.gov/geo/tiger/TIGER2012/COUNTY/tl_2012_us_county.zip
I've not added these files to the working-open-data repo because I don't want large data files in the repo if possible.
import numpy as np
import matplotlib.pyplot as plt
with warnings.catch_warnings():
warnings.simplefilter("ignore")
from mpl_toolkits.basemap import Basemap
# you can figure out where the data directory for basemap is
import os
import inspect
basemap_data_dir = os.path.join(os.path.dirname(inspect.getfile(Basemap)), "data")
print basemap_data_dir, os.path.exists(os.path.join(basemap_data_dir,"UScounties.shp"))
C:\Python27\lib\site-packages\mpl_toolkits\basemap\data True
!ls $basemap_data_dir
'ls' is not recognized as an internal or external command, operable program or batch file.
# https://code.google.com/p/pyshp/
import shapefile
# this is my git clone of https://github.com/matplotlib/basemap --> these files will be in the PiCloud basemap_data_dir
if os.path.exists(os.path.join(basemap_data_dir,"UScounties.shp")):
shpf = shapefile.Reader(os.path.join(basemap_data_dir,"UScounties"))
else:
# put in your path
#shpf = shapefile.Reader("/Users/raymondyee/Dropbox/WwoD13/tl_2012_us_county")
shpf = shapefile.Reader("/Users/raymondyee/C/src/basemap/lib/mpl_toolkits/basemap/data/UScounties")
shapes = shpf.shapes()
records = shpf.records()
for x in records:
print x
break
for y in shapes:
print y
break
['01', '029', '01029', 'AL', 'Cleburne', 'County'] <shapefile._Shape instance at 0x0573C7B0>
shpf.fields
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-10-335ee161db96> in <module>() 1 shpf.fields ----> 2 shapes.parts AttributeError: 'list' object has no attribute 'parts'
zip([f[0] for f in shpf.fields[1:]], records[0])
[('STATE_FIPS', '01'), ('COUNTY_FIP', '029'), ('FIPS', '01029'), ('STATE', 'AL'), ('NAME', 'Cleburne'), ('LSAD', 'County')]
# just CA
from itertools import islice, izip
len([r for r in islice(records,None) if r[0] == '06'])
58
def zip_filter_by_state(records, shapes, included_states=None):
# by default, no filtering
# included_states is a list of states fips prefixes
for (record, state) in izip(records, shapes):
if record[0] in included_states:
yield (record, state)
list(zip_filter_by_state(records, shapes, ['06']))
[(['06', '035', '06035', 'CA', 'Lassen', 'County'], <shapefile._Shape instance at 0x0879FCD8>), (['06', '049', '06049', 'CA', 'Modoc', 'County'], <shapefile._Shape instance at 0x0879FD28>), (['06', '075', '06075', 'CA', 'San Francisco', 'County'], <shapefile._Shape instance at 0x0879FD78>), (['06', '083', '06083', 'CA', 'Santa Barbara', 'County'], <shapefile._Shape instance at 0x0879FDC8>), (['06', '091', '06091', 'CA', 'Sierra', 'County'], <shapefile._Shape instance at 0x0879FE18>), (['06', '095', '06095', 'CA', 'Solano', 'County'], <shapefile._Shape instance at 0x0879FE68>), (['06', '101', '06101', 'CA', 'Sutter', 'County'], <shapefile._Shape instance at 0x0879FEB8>), (['06', '105', '06105', 'CA', 'Trinity', 'County'], <shapefile._Shape instance at 0x0879FF08>), (['06', '111', '06111', 'CA', 'Ventura', 'County'], <shapefile._Shape instance at 0x0879FF58>), (['06', '115', '06115', 'CA', 'Yuba', 'County'], <shapefile._Shape instance at 0x0879FFA8>), (['06', '033', '06033', 'CA', 'Lake', 'County'], <shapefile._Shape instance at 0x0898E030>), (['06', '071', '06071', 'CA', 'San Bernardino', 'County'], <shapefile._Shape instance at 0x0898E080>), (['06', '087', '06087', 'CA', 'Santa Cruz', 'County'], <shapefile._Shape instance at 0x0898E0D0>), (['06', '001', '06001', 'CA', 'Alameda', 'County'], <shapefile._Shape instance at 0x0898E120>), (['06', '003', '06003', 'CA', 'Alpine', 'County'], <shapefile._Shape instance at 0x0898E170>), (['06', '005', '06005', 'CA', 'Amador', 'County'], <shapefile._Shape instance at 0x0898E1C0>), (['06', '007', '06007', 'CA', 'Butte', 'County'], <shapefile._Shape instance at 0x0898E210>), (['06', '009', '06009', 'CA', 'Calaveras', 'County'], <shapefile._Shape instance at 0x0898E260>), (['06', '011', '06011', 'CA', 'Colusa', 'County'], <shapefile._Shape instance at 0x0898E2B0>), (['06', '013', '06013', 'CA', 'Contra Costa', 'County'], <shapefile._Shape instance at 0x0898E300>), (['06', '015', '06015', 'CA', 'Del Norte', 'County'], <shapefile._Shape instance at 0x0898E350>), (['06', '017', '06017', 'CA', 'El Dorado', 'County'], <shapefile._Shape instance at 0x0898E3A0>), (['06', '019', '06019', 'CA', 'Fresno', 'County'], <shapefile._Shape instance at 0x0898E3F0>), (['06', '021', '06021', 'CA', 'Glenn', 'County'], <shapefile._Shape instance at 0x0898E440>), (['06', '023', '06023', 'CA', 'Humboldt', 'County'], <shapefile._Shape instance at 0x0898E490>), (['06', '025', '06025', 'CA', 'Imperial', 'County'], <shapefile._Shape instance at 0x0898E4E0>), (['06', '027', '06027', 'CA', 'Inyo', 'County'], <shapefile._Shape instance at 0x0898E530>), (['06', '029', '06029', 'CA', 'Kern', 'County'], <shapefile._Shape instance at 0x0898E580>), (['06', '031', '06031', 'CA', 'Kings', 'County'], <shapefile._Shape instance at 0x0898E5D0>), (['06', '037', '06037', 'CA', 'Los Angeles', 'County'], <shapefile._Shape instance at 0x0898E620>), (['06', '039', '06039', 'CA', 'Madera', 'County'], <shapefile._Shape instance at 0x0898E670>), (['06', '041', '06041', 'CA', 'Marin', 'County'], <shapefile._Shape instance at 0x0898E6C0>), (['06', '043', '06043', 'CA', 'Mariposa', 'County'], <shapefile._Shape instance at 0x0898E710>), (['06', '045', '06045', 'CA', 'Mendocino', 'County'], <shapefile._Shape instance at 0x0898E760>), (['06', '047', '06047', 'CA', 'Merced', 'County'], <shapefile._Shape instance at 0x0898E7B0>), (['06', '051', '06051', 'CA', 'Mono', 'County'], <shapefile._Shape instance at 0x0898E800>), (['06', '053', '06053', 'CA', 'Monterey', 'County'], <shapefile._Shape instance at 0x0898E850>), (['06', '055', '06055', 'CA', 'Napa', 'County'], <shapefile._Shape instance at 0x0898E8A0>), (['06', '057', '06057', 'CA', 'Nevada', 'County'], <shapefile._Shape instance at 0x0898E8F0>), (['06', '059', '06059', 'CA', 'Orange', 'County'], <shapefile._Shape instance at 0x0898E940>), (['06', '061', '06061', 'CA', 'Placer', 'County'], <shapefile._Shape instance at 0x0898E990>), (['06', '063', '06063', 'CA', 'Plumas', 'County'], <shapefile._Shape instance at 0x0898E9E0>), (['06', '065', '06065', 'CA', 'Riverside', 'County'], <shapefile._Shape instance at 0x0898EA30>), (['06', '067', '06067', 'CA', 'Sacramento', 'County'], <shapefile._Shape instance at 0x0898EA80>), (['06', '069', '06069', 'CA', 'San Benito', 'County'], <shapefile._Shape instance at 0x0898EAD0>), (['06', '073', '06073', 'CA', 'San Diego', 'County'], <shapefile._Shape instance at 0x0898EB20>), (['06', '077', '06077', 'CA', 'San Joaquin', 'County'], <shapefile._Shape instance at 0x0898EB70>), (['06', '079', '06079', 'CA', 'San Luis Obispo', 'County'], <shapefile._Shape instance at 0x0898EBC0>), (['06', '081', '06081', 'CA', 'San Mateo', 'County'], <shapefile._Shape instance at 0x0898EC10>), (['06', '085', '06085', 'CA', 'Santa Clara', 'County'], <shapefile._Shape instance at 0x0898EC60>), (['06', '089', '06089', 'CA', 'Shasta', 'County'], <shapefile._Shape instance at 0x0898ECB0>), (['06', '093', '06093', 'CA', 'Siskiyou', 'County'], <shapefile._Shape instance at 0x0898ED00>), (['06', '097', '06097', 'CA', 'Sonoma', 'County'], <shapefile._Shape instance at 0x0898ED50>), (['06', '099', '06099', 'CA', 'Stanislaus', 'County'], <shapefile._Shape instance at 0x0898EDA0>), (['06', '103', '06103', 'CA', 'Tehama', 'County'], <shapefile._Shape instance at 0x0898EDF0>), (['06', '107', '06107', 'CA', 'Tulare', 'County'], <shapefile._Shape instance at 0x0898EE40>), (['06', '109', '06109', 'CA', 'Tuolumne', 'County'], <shapefile._Shape instance at 0x0898EE90>), (['06', '113', '06113', 'CA', 'Yolo', 'County'], <shapefile._Shape instance at 0x0898EEE0>)]
len(shapes)
3221
# http://www.geophysique.be/2013/02/12/matplotlib-basemap-tutorial-10-shapefiles-unleached-continued/
#
# BaseMap example by geophysique.be
# tutorial 10
import time
import os
import inspect
import numpy as np
import matplotlib.pyplot as plt
from itertools import islice, izip
from mpl_toolkits.basemap import Basemap
### PARAMETERS FOR MATPLOTLIB :
import matplotlib as mpl
mpl.rcParams['font.size'] = 10.
mpl.rcParams['font.family'] = 'Comic Sans MS'
mpl.rcParams['axes.labelsize'] = 8.
mpl.rcParams['xtick.labelsize'] = 6.
mpl.rcParams['ytick.labelsize'] = 6.
fig = plt.figure(figsize=(11.7,8.3))
#Custom adjust of the subplots
plt.subplots_adjust(left=0.05,right=0.95,top=0.90,bottom=0.05,wspace=0.15,hspace=0.05)
ax = plt.subplot(111)
#Let's create a basemap of USA
#x1 = -180.
#x2 = -62.
#y1 = 18.
#y2 = 68.
x1 = -128.
x2 = -111.
y1 = 25.
y2 = 44.5
#m = Basemap(llcrnrlon=-128,llcrnrlat=25,urcrnrlon=-111.0,urcrnrlat=44.5,
# resolution='i',projection='cass',lon_0=-120.36,lat_0=37.7)
m = Basemap(resolution='i',projection='merc', llcrnrlat=y1,urcrnrlat=y2,llcrnrlon=x1,urcrnrlon=x2,lat_ts=(y1+y2)/2)
#m.drawcountries(linewidth=0.5)
#m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(y1,y2,2.),labels=[1,0,0,0],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw parallels
m.drawmeridians(np.arange(x1,x2,2.),labels=[0,0,0,1],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw meridians
def zip_filter_by_state(records, shapes, included_states=None):
# by default, no filtering
# included_states is a list of states fips prefixes
for (record, state) in izip(records, shapes):
if record[0] in included_states:
yield (record, state)
from matplotlib.collections import LineCollection
from matplotlib import cm
import shapefile
basemap_data_dir = os.path.join(os.path.dirname(inspect.getfile(Basemap)), "data")
# this is my git clone of https://github.com/matplotlib/basemap --> these files will be in the PiCloud basemap_data_dir
if os.path.exists(os.path.join(basemap_data_dir,"UScounties.shp")):
shpf = shapefile.Reader(os.path.join(basemap_data_dir,"UScounties"))
else:
# put in your path
#shpf = shapefile.Reader("/Users/raymondyee/Dropbox/WwoD13/tl_2012_us_county")
shpf = shapefile.Reader("/Users/raymondyee/C/src/basemap/lib/mpl_toolkits/basemap/data/UScounties")
shapes = shpf.shapes()
records = shpf.records()
# show only CA and AK (for example)
for record, shape in zip_filter_by_state(records, shapes, ['06']):# removed 02 state , we just need california
lons,lats = zip(*shape.points) # lat long for each county
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
#time.sleep(1) # introduce a delay of 1 second
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
color = 0
if record[2] == '06071': # logic for mapping counties to color gradients - high lighting san-bernandino county
color = 0.1
else:
color = 0.5
lines.set_facecolors(cm.jet(color)) #np.random.rand(1)
lines.set_edgecolors('k')
lines.set_linewidth(0.1)
ax.add_collection(lines) # draw step
time.sleep(1)
plt.savefig('tutorial10.png',dpi=300)
plt.show()
# read deb's dataframe
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
# llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
# are the lat/lon values of the lower left and upper right corners
# of the map.
# resolution = 'i' means use intermediate resolution coastlines.
# lon_0, lat_0 are the central longitude and latitude of the projection.
#m = Basemap(llcrnrlon=-10.5,llcrnrlat=49.5,urcrnrlon=3.5,urcrnrlat=59.5,
# resolution='i',projection='cass',lon_0=-4.36,lat_0=54.7)
m = Basemap(llcrnrlon=-128,llcrnrlat=25,urcrnrlon=-111.0,urcrnrlat=44.5,
resolution='i',projection='cass',lon_0=-120.36,lat_0=37.7)
# can get the identical map this way (by specifying width and
# height instead of lat/lon corners)
#m = Basemap(width=891185,height=1115557,\
# resolution='i',projection='cass',lon_0=-4.36,lat_0=54.7)
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
#m.drawparallels(np.arange(-40,61.,2.))
#m.drawmeridians(np.arange(-20.,21.,2.))
m.drawmapboundary(fill_color='aqua')
plt.title("Cassini Projection")
plt.show()