<jupyter_start><jupyter_text>Example of using ShapefilesI learned this by reading through the webpages linked in the comments throughout thisPlease do reach out if you have any questions about plotting with these!-<NAME><jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
# from: https://github.com/GeospatialPython/pyshp
# install via "pip install pyshp"
import shapefile
# from: https://www.seattle.gov/Documents/Departments/OPCD/Demographics/GeographicFilesandMaps/KingCountyTractsShapefiles.zip
sf = shapefile.Reader('../KingCountyTractsShapefiles/kc_tract_10.shp')

print(str(len(sf.shapes())) + " shapes")

# from: https://chrishavlin.wordpress.com/2016/11/16/shapefiles-tutorial/
plt.figure(figsize=(14,14))
ax = plt.axes() # add the axes
ax.set_aspect('equal')

for shape in list(sf.iterShapes()):
    npoints=len(shape.points) # total points
    nparts = len(shape.parts) # total parts

    if nparts == 1:
        x_lon = np.zeros((len(shape.points),1))
        y_lat = np.zeros((len(shape.points),1))
        for ip in range(len(shape.points)):
            x_lon[ip] = shape.points[ip][0]
            y_lat[ip] = shape.points[ip][1]
        plt.plot(x_lon, y_lat, c='k')

    else: # loop over parts of each shape, plot separately
        for ip in range(nparts): # loop over parts, plot separately
            i0=shape.parts[ip]
            if ip < nparts-1:
                i1 = shape.parts[ip+1]-1
            else:
                i1 = npoints

        seg=shape.points[i0:i1+1]
        x_lon = np.zeros((len(seg),1))
        y_lat = np.zeros((len(seg),1))
        for ip in range(len(seg)):
            x_lon[ip] = seg[ip][0]
            y_lat[ip] = seg[ip][1]

        plt.plot(x_lon, y_lat, c='k')

plt.figure(figsize=(14,14))
ax = plt.axes() # add the axes
ax.set_aspect('equal')

for shape in list(sf.iterShapes()):
    npoints=len(shape.points) # total points
    nparts = len(shape.parts) # total parts

    if nparts == 1:
        x_lon = np.zeros((len(shape.points),1))
        y_lat = np.zeros((len(shape.points),1))
        for ip in range(len(shape.points)):
            x_lon[ip] = shape.points[ip][0]
            y_lat[ip] = shape.points[ip][1]
            
        plt.plot(x_lon, y_lat, c='k')
        plt.fill(x_lon, y_lat, alpha=0.5)

    else: # loop over parts of each shape, plot separately
        for ip in range(nparts): # loop over parts, plot separately
            i0=shape.parts[ip]
            if ip < nparts-1:
                i1 = shape.parts[ip+1]-1
            else:
                i1 = npoints

        seg=shape.points[i0:i1+1]
        x_lon = np.zeros((len(seg),1))
        y_lat = np.zeros((len(seg),1))
        for ip in range(len(seg)):
            x_lon[ip] = seg[ip][0]
            y_lat[ip] = seg[ip][1]

        plt.plot(x_lon, y_lat, c='k')
        plt.fill(x_lon, y_lat, alpha=0.5)<jupyter_output><empty_output><jupyter_text>NOTE:the units for this map are... wtf. They are in "linear feet" from a central point (according to the note on the website I downloaded them from), and I don't quite know what it is.You could find the middle of this map, subtract it off, and convert it to miles or km or something useful!<jupyter_code># let's eyeball a middle
x0 = 1355000
y0 = 152000

ft_to_km = 0.0003048 # accroding to Google

plt.figure(figsize=(14,14))
ax = plt.axes() # add the axes
ax.set_aspect('equal')

for shape in list(sf.iterShapes()):
    npoints=len(shape.points) # total points
    nparts = len(shape.parts) # total parts

    if nparts == 1:
        x_lon = np.zeros((len(shape.points),1))
        y_lat = np.zeros((len(shape.points),1))
        for ip in range(len(shape.points)):
            x_lon[ip] = shape.points[ip][0]
            y_lat[ip] = shape.points[ip][1]
            
        plt.plot((x_lon - x0)*ft_to_km, (y_lat - y0)*ft_to_km, c='k')
        plt.fill((x_lon - x0)*ft_to_km, (y_lat - y0)*ft_to_km, alpha=0.5)

    else: # loop over parts of each shape, plot separately
        for ip in range(nparts): # loop over parts, plot separately
            i0=shape.parts[ip]
            if ip < nparts-1:
                i1 = shape.parts[ip+1]-1
            else:
                i1 = npoints

        seg=shape.points[i0:i1+1]
        x_lon = np.zeros((len(seg),1))
        y_lat = np.zeros((len(seg),1))
        for ip in range(len(seg)):
            x_lon[ip] = seg[ip][0]
            y_lat[ip] = seg[ip][1]

        plt.plot((x_lon - x0)*ft_to_km, (y_lat - y0)*ft_to_km, c='k')
        plt.fill((x_lon - x0)*ft_to_km, (y_lat - y0)*ft_to_km, alpha=0.5)

plt.xlabel('km')
plt.ylabel('km')
plt.title('King County, US 2010 Census Tracts')
# a neat file to connect Census blocks into Neighborhoods!
# from this page: https://www.seattle.gov/opcd/population-and-demographics/geographic-files-and-maps#2010census
# https://www.seattle.gov/Documents/Departments/OPCD/Demographics/GeographicFilesandMaps/SeattleCensusBlocksandNeighborhoodCorrelationFile.xlsx
# from: https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2017&layergroup=Census+Tracts
sf = shapefile.Reader('../tl_2017_53_tract/tl_2017_53_tract.shp')

print(str(len(sf.shapes())) + " shapes")

# from: https://chrishavlin.wordpress.com/2016/11/16/shapefiles-tutorial/
plt.figure(figsize=(14,10))
ax = plt.axes() # add the axes
ax.set_aspect('auto')

for shape in list(sf.iterShapes()):
    npoints=len(shape.points) # total points
    nparts = len(shape.parts) # total parts

    if nparts == 1:
        x_lon = np.zeros((len(shape.points),1))
        y_lat = np.zeros((len(shape.points),1))
        for ip in range(len(shape.points)):
            x_lon[ip] = shape.points[ip][0]
            y_lat[ip] = shape.points[ip][1]
        plt.plot(x_lon, y_lat, c='k', lw=0.5)
        plt.fill(x_lon, y_lat, alpha=0.5)

    else: # loop over parts of each shape, plot separately
        for ip in range(nparts): # loop over parts, plot separately
            i0=shape.parts[ip]
            if ip < nparts-1:
                i1 = shape.parts[ip+1]-1
            else:
                i1 = npoints

        seg=shape.points[i0:i1+1]
        x_lon = np.zeros((len(seg),1))
        y_lat = np.zeros((len(seg),1))
        for ip in range(len(seg)):
            x_lon[ip] = seg[ip][0]
            y_lat[ip] = seg[ip][1]

        plt.plot(x_lon, y_lat, c='k', lw=0.5)
        plt.fill(x_lon, y_lat, alpha=0.5)

plt.xlabel('Lon (deg)')
plt.ylabel('Lat (deg)')
plt.title('WA, US 2010 Census Tracts')<jupyter_output>1458 shapes
{'fraction_non_alphanumeric': 0.1456953642384106, 'fraction_numerical': 0.028492222393346683, 'mean_word_length': 2.6751556310130162, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 9, 'https://': 7, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '436935', 'license': 'MIT', 'n_tokens_mistral': 2505, 'n_tokens_neox': 2355, 'n_words': 506, 'path': 'shapefile_example.ipynb', 'repo_name': 'jradavenport/playing_with_shapefiles', 'size': '6500'}