This notebook goes with the blog post of the same name, published on 10 August 2017.
We're going to load a shapefile containing some well data. Then we'll change its CRS, make a new attribute, and save a new shapefile.
We'll lean on geopandas for help, but we'll also inspect the file with fiona, a lower-level library that geopandas uses under the hood.
Install geopandas and its dependencies (like gdal, proj, and fiona) with
conda install geopandas
conda install fiona # I had to do this too to get fiona to work.
In [1]:
import numpy as np
import fiona
import matplotlib.pyplot as plt
import folium
%matplotlib inline
In [2]:
import pprint
with fiona.open('../data/offshore_wells_2011_Geographic_NAD27.shp') as src:
pprint.pprint(src[0])
In [3]:
import geopandas as gpd
Load our data into a GeoDataFrame (gdf):
In [4]:
gdf = gpd.read_file('../data/offshore_wells_2011_Geographic_NAD27.shp')
In [5]:
gdf.head()
Out[5]:
In [6]:
gdf.plot()
Out[6]:
Let's look at the first 5 rows of the geometry attribute.
In [7]:
gdf.geometry[:5]
Out[7]:
Notice we're in lat, lon:
In [8]:
print(gdf.crs)
Visiting EPSG 4267 tells us the datum is NAD27.
Let's cast the SHP to a new CRS: EPSG 26920:
In [9]:
gdf = gdf.to_crs({'init': 'epsg:26920'})
Now we're in a UTM coordinate system: UTM Zone 20N, with a NAD83 datum:
In [10]:
gdf.geometry[:5]
Out[10]:
Now we'll also add a new attribute with the two-way seismic travel time to the seafloor (in milliseconds).
In [11]:
gdf['seafl_twt'] = 2 * 1000 * gdf.Water_Dept / 1485
In [12]:
gdf.head()
Out[12]:
We can also get a statistical summary of the data frame:
In [13]:
gdf.describe()
Out[13]:
In [14]:
gdf.to_file('../data/offshore_wells_2011_UTM20_NAD83.shp')
In [15]:
ls ../data/*.shp
In [16]:
# Must be geographic coords, so casting to WGS84.
gdf = gdf.to_crs({'init': 'epsg:4326'})
with open('../data/offshore_wells.geojson', 'w') as f:
f.write(gdf.to_json())
In [17]:
import folium
# Must be geographic coords, so casting to WGS84.
gdf = gdf.to_crs({'init': 'epsg:4326'})
# Make the map, add the features via GeoJSON.
mymap = folium.Map(location=[45, -62], zoom_start=7)
features = folium.features.GeoJson(gdf.to_json())
mymap.add_child(features)
mymap
Out[17]:
© Agile Geoscience 2016