Map of New Mexico's Tribal lands and grocery stores

For this map I used data collected using ArcGIS, specifically data regarding New Mexico's grcoery stores and New Mexico's tribal lands.

In [33]:
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
In [ ]:
 
In [34]:
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'
In [35]:
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.

In [36]:
glaciers.head()
Out[36]:
BINOM_BEG C CM COUNTY C_rt C_se C_si C_uc EST_NAME FID_ ... R_rm_c R_rm_i R_rm_j R_ur SOURCE U UID_JOIN_F UID_JOIN_L WD geometry
0 999999.0 0.0 0.0 31.0 0.0 0.0 0.0 0.0 Albertsons - No902 0 ... None None None None D5 0.0 5538.0 5538.0 0.0 POINT (-108.71737 35.529602)
1 999999.0 0.0 0.0 27.0 0.0 0.0 0.0 0.0 Wal Mart Super Center - No851 0 ... None None None None D1 0.0 788.0 788.0 0.0 POINT (-105.547561 33.35897400000005)
2 999999.0 0.0 0.0 9.0 0.0 0.0 0.0 0.0 Albertsons - No916 0 ... None None None None D4 0.0 4077.0 4077.0 0.0 POINT (-103.196493 34.41744100000002)
3 999999.0 0.0 0.0 17.0 0.0 0.0 0.0 0.0 Albertsons - No 914 0 ... None None None None D3 0.0 3580.0 3580.0 0.0 POINT (-108.255132 32.78739199999997)
4 999999.0 1.0 0.0 49.0 1.0 0.0 0.0 0.0 Sunflower Farmers Markets 0 ... None None None None D2 0.0 2023.0 2023.0 0.0 POINT (-105.951133 35.69142899999996)

5 rows × 69 columns

Here is the head of the tribal GeoPandas DataFrame.

In [37]:
tribal.head()
Out[37]:
AIANNHCE00 AIANNHID00 AIANNHR00 CLASSFP00 COMPTYP00 FID_ FUNCSTAT00 LSAD00 MTFCC00 NAME00 NAMELSAD00 geometry
0 2430 2430R F D8 R 0 A 86 G2101 Navajo Nation Navajo Nation Reservation (POLYGON ((-107.7396288784021 36.6471030979193...
1 4470 4470R F D8 R 0 A 86 G2101 Ute Mountain Ute Mountain Reservation POLYGON ((-108.4628379999998 36.85005099999989...
2 4785 4785R F D8 R 0 A 86 G2101 Zuni Zuni Reservation (POLYGON ((-109.0460690081911 34.8578950239504...
3 0010 0010R F D8 R 0 A 84 G2101 Acoma Acoma Pueblo (POLYGON ((-107.5585460623441 35.0739959704343...
4 2400 2400R F D8 R 0 A 84 G2101 Nambe Nambe Pueblo POLYGON ((-105.9289460000002 35.82681299999997...

Now, here is the head of the counties GeoPandas dataframe

In [38]:
counties.head()
Out[38]:
ALAND AWATER CBSAFP CLASSFP COUNTYFP COUNTYNS CSAFP FID_ FUNCSTAT GEOID INTPTLAT INTPTLON LSAD METDIVFP MTFCC NAME NAMELSAD STATEFP geometry
0 7.679550e+09 431978.0 19700 H1 029 00933057 None 0 A 35029 +32.1844822 -107.7466387 06 None G4020 Luna Luna County 35 POLYGON ((-108.2297870000003 32.20716300000004...
1 8.901430e+09 23652092.0 None H1 023 00929106 None 0 A 35023 +31.8996575 -108.7457292 06 None G4020 Hidalgo Hidalgo County 35 POLYGON ((-109.0486389999999 31.92091300000003...
2 1.026056e+10 15349508.0 43500 H1 017 00915980 None 0 A 35017 +32.7320870 -108.3815043 06 None G4020 Grant Grant County 35 POLYGON ((-108.6594840000004 33.2093790000001,...
3 3.006531e+09 16471227.0 10740 H1 001 01702363 None 0 A 35001 +35.0540019 -106.6690645 06 None G4020 Bernalillo Bernalillo County 35 POLYGON ((-106.3731929999999 34.87033799999998...
4 1.517950e+10 91459429.0 21580 H1 039 01702368 492 0 A 35039 +36.5096687 -106.6939829 06 None G4020 Rio Arriba Rio Arriba County 35 POLYGON ((-106.9855130000004 37.00012999999998...

Population density of all the tribes from the tribal land GeoDataframe

In [39]:
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.

In [40]:
y = glaciers.LAT.values
x = glaciers.LONG_.values

Composing the basemap of New Mexico:

In [41]:
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.

In [42]:
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()
['CLASSFP00', 'FID_', 'MTFCC00', 'LSAD00', 'NAME00', 'AIANNHR00', 'RINGNUM', 'AIANNHID00', 'FUNCSTAT00', 'NAMELSAD00', 'COMPTYP00', 'AIANNHCE00', 'SHAPENUM']

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.

In [43]:
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()
['CLASSFP00', 'FID_', 'MTFCC00', 'LSAD00', 'NAME00', 'AIANNHR00', 'RINGNUM', 'AIANNHID00', 'FUNCSTAT00', 'NAMELSAD00', 'COMPTYP00', 'AIANNHCE00', 'SHAPENUM']

As we can see, the majority of grocery stores in New Mexico are not located within the tribal lands.

In [ ]: