{"id":224,"date":"2014-06-17T15:09:24","date_gmt":"2014-06-17T15:09:24","guid":{"rendered":"http:\/\/blogs.law.harvard.edu\/rprasad\/?p=224"},"modified":"2014-06-17T15:53:46","modified_gmt":"2014-06-17T15:53:46","slug":"shapefile-reprojected-shp-geojson-leaflet","status":"publish","type":"post","link":"https:\/\/archive.blogs.harvard.edu\/rprasad\/2014\/06\/17\/shapefile-reprojected-shp-geojson-leaflet\/","title":{"rendered":"Shapefile -> Reprojected shp -> GeoJson -> Leaflet"},"content":{"rendered":"

\"Mozilla_Firefox\"<\/a><\/p>\n

Here’s a not-too-pretty proof of concept script that hastily glues together several blog posts, as noted in the code.<\/p>\n

The script takes the following steps:<\/p>\n

    \n
  1. Takes the path to a .shp file<\/li>\n
  2. Reprojects the .shp file to 4326<\/a> (without checking if it’s already in 4326–DOH!) (via Python GDAL\/OGR Cookbook<\/a>)<\/li>\n
  3. Converts the .shp to GEOJSON<\/a> (by Shapefile to GeoJSON<\/a>)<\/li>\n
  4. Creates a Leaflet page using Folium (using Folium<\/a>)<\/li>\n
  5. Does all this using way too many libraries:) <\/li>\n<\/ol>\n

    from __future__ import print_function\nimport os, sys\nimport json\n\nfrom osgeo import ogr, osr\nimport shapefile\nimport folium\n\ndef msg(s): print (s)\ndef dashes(): msg(40*'-')\ndef msgt(s): dashes(); msg(s); dashes()\ndef msgx(s): dashes(); msg('ERROR'); msg(s); dashes(); sys.exit(0)\n\n\ndef get_output_fname(fname, new_suffix):\n    fparts = fname.split('.')\n    if len(fparts[-1]) == 3:\n        return '.'.join(fparts[:-1]) + new_suffix + '.' + fparts[-1]\n\n    return fname + new_suffix\n\ndef reproject_to_4326(shape_fname):\n    \"\"\"Re-project the shapefile to a 4326. From the Python GDAL\/OGR Cookbook\n\n    Source: http:\/\/pcjericks.github.io\/py-gdalogr-cookbook\/projection.html#reproject-a-layer\n\n    :param shape_fname: full file path to a shapefile (.shp)\n    :returns: full file path to a shapefile reprojected as 4326\n    \"\"\"\n    if not os.path.isfile(shape_fname):\n        msgx('File not found: %s' % shape_fname)\n\n    driver = ogr.GetDriverByName('ESRI Shapefile')\n    inDataSet = driver.Open(shape_fname)\n\n    # input SpatialReference\n    inLayer = inDataSet.GetLayer()\n    inSpatialRef = inLayer.GetSpatialRef()\n\n    # output SpatialReference\n    outSpatialRef = osr.SpatialReference()\n    outSpatialRef.ImportFromEPSG(4326)\n\n    # create the CoordinateTransformation\n    coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)\n\n    # create the output layer\n    outputShapefile = get_output_fname(shape_fname, '_4326')\n    #msg('output file: %s' % outputShapefile)\n\n    if os.path.exists(outputShapefile):\n        driver.DeleteDataSource(outputShapefile)\n    outDataSet = driver.CreateDataSource(outputShapefile)\n    outLayer = outDataSet.CreateLayer(\"basemap_4326\", geom_type=ogr.wkbMultiPolygon)\n\n    # add fields\n    inLayerDefn = inLayer.GetLayerDefn()\n    for i in range(0, inLayerDefn.GetFieldCount()):\n        fieldDefn = inLayerDefn.GetFieldDefn(i)\n        outLayer.CreateField(fieldDefn)\n\n    # get the output layer's feature definition\n    outLayerDefn = outLayer.GetLayerDefn()\n\n    # loop through the input features\n    inFeature = inLayer.GetNextFeature()\n    while inFeature:\n        # get the input geometry\n        geom = inFeature.GetGeometryRef()\n        # reproject the geometry\n        geom.Transform(coordTrans)\n        # create a new feature\n        outFeature = ogr.Feature(outLayerDefn)\n        # set the geometry and attribute\n        outFeature.SetGeometry(geom)\n        for i in range(0, outLayerDefn.GetFieldCount()):\n            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))\n        # add the feature to the shapefile\n        outLayer.CreateFeature(outFeature)\n        # destroy the features and get the next input feature\n        outFeature.Destroy()\n        inFeature.Destroy()\n        inFeature = inLayer.GetNextFeature()\n\n    # close the shapefiles\n    inDataSet.Destroy()\n    outDataSet.Destroy()\n\n\n    msg('output file: %s' % outputShapefile)\n    return outputShapefile\n\ndef convert_shp_to_geojson(shape_fname):\n    \"\"\"Using the pyshp library, https:\/\/github.com\/GeospatialPython\/pyshp, convert the shapefile to JSON\n\n    Code is from this example:  http:\/\/geospatialpython.com\/2013\/07\/shapefile-to-geojson.html\n\n    :param shape_fname: full file path to a shapefile (.shp)\n    :returns: full file path to a GEOJSON representation of the shapefile\n\n    (recheck\/redo using gdal)\n    \"\"\"\n    if not os.path.isfile(shape_fname):\n        msgx('File not found: %s' % shape_fname)\n\n    # Read the shapefile\n    try: \n        reader = shapefile.Reader(shape_fname)\n    except:\n        msgx('Failed to read shapefile: %s' % shape_fname)\n\n    fields = reader.fields[1:]\n    field_names = [field[0] for field in fields]\n    output_buffer = []\n    for sr in reader.shapeRecords():\n          atr = dict(zip(field_names, sr.record))\n          geom = sr.shape.__geo_interface__\n          output_buffer.append(dict(type=\"Feature\", geometry=geom, properties=atr))\n\n    # write the GeoJSON file\n    out_fname = os.path.join('page-out', os.path.basename(shape_fname).replace('.shp', '.json'))\n\n    geojson = open(out_fname, \"w\")\n    geojson.write(json.dumps({\"type\": \"FeatureCollection\",\"features\": output_buffer}, indent=2) + \"n\")\n    geojson.close()    \n    msg('file written: %s' % out_fname)\n\n    return out_fname\n\n\ndef make_leaflet_page(geojson_file, ouput_html_fname):\n    \"\"\"Using folium, make an HTML page using GEOJSON input \n        examples: https:\/\/github.com\/wrobstory\/folium\n\n    :param geojson_file: full file path to a GEO JSON file\n    :param ouput_html_fname: name of HTML file to create (will only use the basename)\n    \"\"\"\n    if not os.path.isfile(geojson_file):\n        msgx('File not found: %s' % geojson_file)\n\n    # Boston \n    mboston = folium.Map(location=[ 42.3267154, -71.1512353])\n    mboston.geo_json(geojson_file)\n\n    #mboston.geo_json('income.json')\n    ouput_html_fname = os.path.basename(ouput_html_fname)\n    mboston.create_map(path=ouput_html_fname)\n    print ('file written', ouput_html_fname)\n\n\nif __name__=='__main__':\n    reprojected_fname = reproject_to_4326('data\/social_disorder_in_boston\/social_disorder_in_boston_yqh.shp')\n    geojson_fname = convert_shp_to_geojson(reprojected_fname)\n    make_leaflet_page(geojson_fname, 'disorder.html')\n<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"
    Here’s a not-too-pretty proof of concept script that hastily glues together several blog posts, as noted in the code. The script takes the following steps: Takes the path to a .shp file Reprojects the .shp file to 4326 (without checking if it’s already in 4326–DOH!) (via Python GDAL\/OGR Cookbook) Converts the .shp to GEOJSON (by … Continue reading Shapefile -> Reprojected shp -> GeoJson -> Leaflet<\/span><\/a><\/p>\n","protected":false},"author":3875,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_feature_clip_id":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2},"jetpack_post_was_ever_published":false},"categories":[1],"tags":[],"class_list":["post-224","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"jetpack_shortlink":"https:\/\/wp.me\/p4JC3p-3C","_links":{"self":[{"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/posts\/224","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/users\/3875"}],"replies":[{"embeddable":true,"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/comments?post=224"}],"version-history":[{"count":12,"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/posts\/224\/revisions"}],"predecessor-version":[{"id":237,"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/posts\/224\/revisions\/237"}],"wp:attachment":[{"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/media?parent=224"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/categories?post=224"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/archive.blogs.harvard.edu\/rprasad\/wp-json\/wp\/v2\/tags?post=224"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}