In this notebook we will:

  • open a Zip file in-memory
  • extract each of the underlying shapefiles
  • convert them to Fiona objects
  • and then plot them.

Extracting a Shapefile from an in-memory folder

Let's try Fiona's virtual file system feature as described here.

In [1]:
import fiona
import os
from zipfile import ZipFile

assets = '/'.join(os.getcwd().split('/')[:-1] + ['assets'])
zipfile = assets + '/shapefiles.zip'
In [3]:
with fiona.drivers():

    for i, layername in enumerate(
            fiona.listlayers(
                '/shapefiles',
                vfs='zip://'+zipfile)):
        with fiona.open(
                '/shapefiles',
                vfs='zip://'+zipfile,
                layer=i) as c:
            print(i, layername, len(c))
(0, u'RoadAsphalt', 2)
(1, u'DitchLine', 40)
(2, u'WaterPipe', 33)
(3, u'boundary', 6)
(4, u'TreeLine', 6)
(5, u'NaturalGas', 53)
(6, u'RoadConcrete', 4)
(7, u'wetlands1', 26258)
(8, u'RoadDriveway', 169)
(9, u'PowerLine', 6)
(10, u'ProposedAccessRoads', 47)
(11, u'Fence', 164)
(12, u'turbines', 49)
(13, u'wetlands2', 32423)
(14, u'RoadGravel', 18)

Beautiful. Now, let's fetch the 'Fence' layer and plot it using Folium into a Leaflet map. First, let's see what is the CRS of our layers.

In [4]:
#Check CRS 
with fiona.drivers():
    with fiona.open(
        '/shapefiles',
        vfs = 'zip://' + zipfile,
        layer = 'Fence') as layer:
        print 'CRS:', layer.crs
CRS: {u'lon_0': -93.5, u'datum': u'NAD83', u'y_0': 0, u'no_defs': True, u'proj': u'lcc', u'x_0': 500000.0000000002, u'units': u'us-ft', u'lat_2': 41.78333333333333, u'lat_1': 40.61666666666667, u'lat_0': 40}

Now, we want to convert this to WGS84 coordinate system to make it easier for Folium to plot. I'll take some code from the WindOps repository and adapt it to this specific case.

In [5]:
import pyproj

def pointTrans(crs,inverse=False):
    proj = pyproj.Proj(crs,preserve_units=True)
    def trans(tpl):
        return proj(tpl[0],tpl[1],inverse=inverse)
    return trans

def customTransform(crs,original,toLatLon=False):
    '''
    Transorms from custom projection to lat,lon
    and viceversa
    crs is a dictionary, original is a geojson-like dict,
    toLatLon: if true will transform to LatLon, otherwise from
    '''
    geometry = original.copy()
    trans = pointTrans(crs,inverse=toLatLon)
    if geometry['type'] == 'Point':
        geometry['coordinates'] = trans(geometry['coordinates'])
    elif geometry['type'] == 'MultiPoint' or geometry['type'] == 'LineString':
        geometry['coordinates'] = map(trans,geometry['coordinates'])
    elif geometry['type'] == 'Polygon' or geometry['type'] == 'MultiLineString':
        elements = []
        for element in geometry['coordinates']:
            elements.append(map(trans,element))
        geometry['coordinates'] = elements
    elif geometry['type'] == 'MultiPolygon':
        elements = []
        for element in geometry['coordinates']:
            subelements = []
            for subelement in element:
                subelements.append(map(trans,subelement))
            elements.append(map(trans, element))
        geometry['coordinates'] = elements
    else:
        raise NameError("Geometry type not supported: " + geometry['type'])
        return
    return geometry


def get_geojson_from_shapefile(layername,zipfile):
    features = []
    with fiona.drivers():
        with fiona.open(
        '/shapefiles',
        vfs = 'zip://' + zipfile,
        layer = layername) as shp:
            crs = shp.crs
            proj = pyproj.Proj(crs)
            for feat in shp:
                feature = feat.copy()
                if not proj.is_latlong():
                    feature['geometry'] = customTransform(shp.crs,feature['geometry'],toLatLon=True)
                features.append(feature)
        geojson = {"type": "FeatureCollection","features": features}
    return geojson
In [6]:
fences = get_geojson_from_shapefile('Fence',zipfile)

We now have the 'Fence' as a GeoJson-like Python dictionary in WGS84 projection. Let's add it to a folium map.

In [7]:
import folium
import json

map_terrain = folium.Map(location=[41.12,-94.11], zoom_start=11)
In [8]:
map_terrain.geo_json(geo_str = json.dumps(fences))
In [9]:
map_terrain
Out[9]:

Excellent, we can see the fences around Macksburg.