For this map I used data collected using ArcGIS, specifically data regarding New Mexico's grcoery stores and New Mexico's tribal lands.
import os
os.chdir('C:\\Users\\jenat\\Documents\\new')
from mpl_toolkits.basemap import Basemap
import matplotlib as mpl
import numpy as np
from descartes import PolygonPatch
import geopandas as gp
import shapefile
from matplotlib.collections import LineCollection
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import shapely
from matplotlib.collections import PatchCollection
import fiona
import shapely.geometry as geometry
import pylab as pl
from __future__ import print_function
from matplotlib.colors import rgb2hex
from matplotlib.patches import Polygon
import os
os.chdir('C:\\Users\\jenat\\Documents\\new')
from mpl_toolkits.basemap import Basemap
import matplotlib as mpl
import numpy as np
from descartes import PolygonPatch
import geopandas as gp
import shapefile
from matplotlib.collections import LineCollection
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from descartes import PolygonPatch
import shapely
earth = 'cornsilk'
from matplotlib.collections import PatchCollection
import fiona
import shapely.geometry as geometry
import pylab as pl
earth = 'cornsilk'
glaciers = gp.GeoDataFrame.from_file('nm_gro - Copy.shp')
counties = gp.GeoDataFrame.from_file('nm_counties - Copy.shp')
tribal = gp.GeoDataFrame.from_file('nm_tl - Copy.shp')
Here is the head of the GeoPandas dataframe regarding New Mexico's grocery stores. I believe this data does NOT include gas stations.
glaciers.head()
Here is the head of the tribal GeoPandas DataFrame.
tribal.head()
Now, here is the head of the counties GeoPandas dataframe
counties.head()
Population density of all the tribes from the tribal land GeoDataframe
popdensity = {
'Navajo Nation Reservation': 138.00,
'Ute Mountain Reservation': 337.35,
'Zuni Reservation': 312.68,
'Acoma Pueblo': 271.40,
'Nambe Pueblo': 209.23,
'Pojoaque Pueblo': 155.18,
'Tesuque Pueblo': 154.87,
'Mescalero Reservation': 114.43,
'Jicarilla Apache Nation Reservation': 107.05,
'San Felipe Pueblo': 105.80,
'San Ildefonso Pueblo': 86.27,
'Santo Domingo Pueblo': 83.85,
'Pueblo de Cochiti': 72.83,
'Laguna Pueblo': 69.03,
'Isleta Pueblo': 67.55,
'Sandia Pueblo': 65.46,
'Santa Clara Pueblo': 63.80,
'Zia Pueblo': 54.59,
'Jemez Pueblo': 53.29,
'Santa Ana Pueblo': 53.20,
'San Juan Pueblo': 51.45,
'Taos Pueblo': 39.61,
'Picuris Pueblo': 39.28,
'Navajo Nation Reservation': 138.00,
'Acoma Pueblo': 271.40,
'Zuni Reservation': 312.68,
'Laguna Pueblo': 69.03,
'Tesuque Pueblo': 154.87,
'Nambe Pueblo': 209.23,
'Zia Pueblo': 54.59,
'Taos Pueblo': 39.61,
'San Felipe/Santa Ana joint-use area': 38.13,
'San Felipe/Santo Domingo joint-use area': 34.20 }
Using the grocery dataframe, I will set my x
and y
values to the latitude
and longitude
values corresponding to the grocery stores locations. I will be using columns from this dataframe which are specificially longitude and latitude.
y = glaciers.LAT.values
x = glaciers.LONG_.values
Composing the basemap of New Mexico:
fig, ax1 = plt.subplots(figsize=(18, 45))
map1 = Basemap(resolution='i',projection='merc',llcrnrlat=31.547, urcrnrlat=37.149, llcrnrlon=-109.0, urcrnrlon=-103.0)
map1.drawmapboundary(fill_color = 'lightskyblue')
map1.fillcontinents(color=earth)
map1.drawcounties(color='purple')
map1.readshapefile('nm_tl - Copy', name='states', drawbounds=True)
plt.title("Map of New Mexico's Tribal lands", fontsize = 20)
plt.show()
This is a map of the state of New Mexico, as well as a shapefile
of the tribal lands which are outlined in black
.
Now, I will colour in all the tribal land polygons.
fig, ax1 = plt.subplots(figsize=(18, 45))
map1 = Basemap(resolution='i',projection='merc',llcrnrlat=31.547, urcrnrlat=37.149, llcrnrlon=-109.0, urcrnrlon=-103.0)
map1.drawmapboundary(fill_color = 'lightskyblue')
map1.fillcontinents(color=earth)
map1.drawcounties(color='purple')
map1.readshapefile('nm_tl - Copy', name='states', drawbounds=True)
###
#colouring in the tribal lands
colors={}
statenames=[]
cmap = plt.cm.hot # use 'hot' colormap
vmin = 0; vmax = 450 # set range.
print(map1.states_info[0].keys())
for shapedict in map1.states_info:
statename = shapedict['NAMELSAD00']
if statename not in ['Cherokee','Sioux']:
pop = popdensity[statename]
# calling colormap with value between 0 and 1 returns
# rgba value. Invert color range (hot colors are high
# population), take sqrt root to spread out colors more.
colors[statename] = cmap(1.-np.sqrt((pop-vmin)/(vmax-vmin)))[:3]
statenames.append(statename)
# cycle through state names, color each one.
ax = plt.gca() # get current axes instance
for nshape,seg in enumerate(map1.states):
if statenames[nshape] not in ['Cherokee','Sioux']:
color = rgb2hex(colors[statenames[nshape]])
poly = Polygon(seg,facecolor=color,edgecolor=color)
ax.add_patch(poly)
plt.title("Map of New Mexico's Tribal lands", fontsize = 30)
plt.show()
Next, I will add a scatter plot, which will indicate in blue the grocery stores located within New Mexico. It should be noted that the grocery store data does not include gas stations.
fig, ax1 = plt.subplots(figsize=(18, 45))
map1 = Basemap(resolution='i',projection='merc',llcrnrlat=31.547, urcrnrlat=37.149, llcrnrlon=-109.0, urcrnrlon=-103.0)
map1.drawmapboundary(fill_color = 'lightskyblue')
map1.fillcontinents(color=earth)
map1.drawcounties(color='purple')
x1,y1 = map1(x, y)
map1.scatter(x1, y1, c = 'b',marker = "o", alpha = 0.9, zorder=10)
map1.readshapefile('nm_tl - Copy', name='states', drawbounds=True)
###
#colouring in the tribal lands
colors={}
statenames=[]
cmap = plt.cm.hot # use 'hot' colormap
vmin = 0; vmax = 450 # set range.
print(map1.states_info[0].keys())
for shapedict in map1.states_info:
statename = shapedict['NAMELSAD00']
if statename not in ['Cherokee','Sioux']:
pop = popdensity[statename]
# calling colormap with value between 0 and 1 returns
# rgba value. Invert color range (hot colors are high
# population), take sqrt root to spread out colors more.
colors[statename] = cmap(1.-np.sqrt((pop-vmin)/(vmax-vmin)))[:3]
statenames.append(statename)
# cycle through state names, color each one.
ax = plt.gca() # get current axes instance
for nshape,seg in enumerate(map1.states):
if statenames[nshape] not in ['Cherokee','Sioux']:
color = rgb2hex(colors[statenames[nshape]])
poly = Polygon(seg,facecolor=color,edgecolor=color)
ax.add_patch(poly)
plt.title("Map of New Mexico's Tribal lands and the amount of grocery store", fontsize = 30)
plt.show()
As we can see, the majority of grocery stores in New Mexico are not located within the tribal lands.