from swmmio.defs.config import ROOT_DIR
import json
import pandas as pd
from geojson import Point, LineString, Polygon, FeatureCollection, Feature
import os
import shutil
[docs]
def change_crs(series, in_crs, to_crs):
"""
Change the projection of a series of coordinates.
Parameters
----------
series : pd.Series, pd.DataFrame, list, or tuple
Series of coordinates to be reprojected.
in_crs : str
The current coordinate reference system of the series.
to_crs : str
The target coordinate reference system to reproject the series to.
Returns
-------
pd.Series, pd.DataFrame, list, or tuple
Object with reprojected coordinates.
Examples
--------
>>> from swmmio.examples import jersey
>>> proj4_str = '+proj=tmerc +lat_0=36.16666666666666 +lon_0=-94.5 +k=0.9999411764705882 +x_0=850000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs'
>>> jersey.crs = proj4_str
>>> change_crs(jersey.nodes.dataframe['coords'], proj4_str, "EPSG:4326")
Name
J3 [(39.236286854940964, -94.64346373821752)]
1 [(39.23851590020802, -94.64756446847099)]
2 [(39.2382157223383, -94.64468629488778)]
3 [(39.23878251491925, -94.64640342340165)]
4 [(39.238353081411915, -94.64603818939938)]
5 [(39.23797714290924, -94.64589184224722)]
J2 [(39.23702605103406, -94.64543916929885)]
J4 [(39.23633648359375, -94.64190240294558)]
J1 [(39.23723558954326, -94.64583338271147)]
Name: coords, dtype: object
"""
try:
import pyproj
from pyproj import Transformer
except ImportError:
raise ImportError('pyproj module needed. get this package here: ',
'https://pypi.python.org/pypi/pyproj')
# SET UP THE TO AND FROM COORDINATE PROJECTION
transformer = Transformer.from_crs(in_crs, to_crs, always_xy=True)
# convert coords in coordinates, vertices, and polygons inp sections
# transform to the typical 'WGS84' coord system
def get_xys(xy_row):
# need to reverse to lat/long after conversion
return [transformer.transform(x, y) for x, y in xy_row]
if isinstance(series, pd.Series):
# unpack the nested coords
xs = [xy[0][0] for xy in series.to_list()]
ys = [xy[0][1] for xy in series.to_list()]
# transform the whole series at once
xs_trans, ys_trans = transformer.transform(xs, ys)
# increase nest level and return pd.Series
return pd.Series(index=series.index, data=[[coords] for coords in zip(xs_trans, ys_trans)])
if isinstance(series, pd.DataFrame):
# zipped_coords = list(zip(series.X, series.Y))
xs_trans, ys_trans = transformer.transform(series.X, series.Y)
df = pd.DataFrame(data=zip(xs_trans, ys_trans), columns=["X", "Y"], index=series.index)
return df
elif isinstance(series, (list, tuple)):
if isinstance(series[0], (list, tuple)):
# unpack the nested coords
xs = [x for x, y in series]
ys = [y for x, y in series]
return list(zip(transformer.transform(xs, ys)))
else:
return transformer.transform(*series)
[docs]
def coords_series_to_geometry(coords, geomtype='linestring', dtype='geojson'):
"""
Convert a series of coordinates to a series of geometry objects.
Parameters
----------
coords : pd.Series
Series of lists of xy coordinates.
geomtype : str, optional
Geometry type of target ('linestring', 'point', 'polygon'). Default is 'linestring'.
dtype : str, optional
Format of geometry objects to be created ('geojson', 'shapely'). Default is 'geojson'.
Returns
-------
pd.Series
Series of geometry objects.
Examples
--------
>>> import swmmio
>>> model = swmmio.Model(MODEL_FULL_FEATURES_XY)
>>> nodes = model.nodes.dataframe
>>> geoms = coords_series_to_geometry(nodes['coords'], geomtype='point')
>>> geoms.iloc[0]
{"coordinates": [2748073.3060000003, 1117746.087], "type": "Point"}
"""
# detect whether LineString or Point should be used
geomtype = geomtype.lower()
if geomtype == 'linestring':
geoms = [LineString(latlngs) for latlngs in coords]
elif geomtype == 'point':
geoms = [Point(latlngs[0]) for latlngs in coords]
elif geomtype == 'polygon':
geoms = [Polygon([latlngs]) for latlngs in coords]
if dtype.lower() == 'shape':
# convert to shapely objects
try:
from shapely.geometry import shape
except ImportError:
raise ImportError('shapely module needed. Install it via GeoPandas with conda: ',
'conda install geopandas')
geoms = [shape(g) for g in geoms]
return pd.Series(index=coords.index, name='geometry', data=geoms)
[docs]
def write_geojson(df, filename=None, geomtype='linestring', drop_na=True):
"""
Convert dataframe with coords series to geojson format.
Parameters
----------
df : pd.DataFrame
Target dataframe.
filename : str, optional
Path of new file to contain geojson.
geomtype : str, optional
Geometry type ('linestring', 'point', 'polygon'). Default is 'linestring'.
drop_na : bool, optional
Whether to remove properties with None values. Default is True.
Returns
-------
geojson.FeatureCollection
GeoJSON FeatureCollection.
Examples
--------
>>> from swmmio.examples import philly
>>> geoj = write_geojson(philly.links.dataframe, drop_na=True)
>>> print(json.dumps(geoj['features'][0]['properties'], indent=2))
{
"InletNode": "J1-025",
"OutletNode": "J1-026",
"Length": 309.456216,
"Roughness": 0.014,
"InOffset": 0,
"OutOffset": 0.0,
"InitFlow": 0,
"MaxFlow": 0,
"Shape": "CIRCULAR",
"Geom1": 1.25,
"Geom2": 0,
"Geom3": 0,
"Geom4": 0,
"Barrels": 1,
"Name": "J1-025.1"
}
>>> print(json.dumps(geoj['features'][0]['geometry'], indent=2))
{
"type": "LineString",
"coordinates": [
[
2746229.223,
1118867.764
],
[
2746461.473,
1118663.257
]
]
}
"""
# CONVERT THE DF INTO JSON
df['Name'] = df.index # add a name column (we wont have the index)
records = json.loads(df.to_json(orient='records'))
# ITERATE THROUGH THE RECORDS AND CREATE GEOJSON OBJECTS
features = []
for rec in records:
coordinates = rec['coords']
del rec['coords'] # delete the coords so they aren't in the properties
if drop_na:
rec = {k: v for k, v in rec.items() if v is not None}
latlngs = coordinates
if geomtype == 'linestring':
geometry = LineString(latlngs)
elif geomtype == 'point':
geometry = Point(latlngs[0])
elif geomtype == 'polygon':
geometry = Polygon([latlngs])
feature = Feature(geometry=geometry, properties=rec)
features.append(feature)
if filename is not None:
with open(filename, 'wb') as f:
f.write(json.dumps(FeatureCollection(features)))
return filename
else:
return FeatureCollection(features)
[docs]
def write_shapefile(df, filename, geomtype='line', prj=None):
"""
create a shapefile given a pandas Dataframe that has coordinate data in a
column called 'coords'.
"""
import shapefile
df['Name'] = df.index
# create a shp file writer object of geom type 'point'
if geomtype == 'point':
w = shapefile.Writer(filename, shapefile.POINT, autoBalance=True)
elif geomtype == 'line':
w = shapefile.Writer(filename, shapefile.POLYLINE, autoBalance=True)
elif geomtype == 'polygon':
w = shapefile.Writer(filename, shapefile.POLYGON, autoBalance=True)
# add the fields
for fieldname in df.columns:
w.field(fieldname, "C")
for k, row in df.iterrows():
w.record(*row.tolist())
if geomtype == 'line':
w.line([row.coords])
if geomtype == 'point':
w.point(*row.coords[0])
if geomtype == 'polygon':
w.poly([row.coords])
w.close()
# add projection data to the shapefile,
if prj is None:
# if not sepcified, the default, projection is used (PA StatePlane)
prj = os.path.join(ROOT_DIR, 'swmmio/defs/default.prj')
prj_filepath = os.path.splitext(filename)[0] + '.prj'
shutil.copy(prj, prj_filepath)
[docs]
def read_shapefile(shp_path):
"""
Read a shapefile into a Pandas dataframe with a 'coords' column holding
the geometry information. This uses the pyshp package
"""
import shapefile
# read file, parse out the records and shapes
sf = shapefile.Reader(shp_path)
fields = [x[0] for x in sf.fields][1:]
records = sf.records()
shps = [s.points for s in sf.shapes()]
# write into a dataframe
df = pd.DataFrame(columns=fields, data=records)
df = df.assign(coords=shps)
return df
[docs]
def centroid_and_bbox_from_coords(coords):
"""
return a bounding box that encapsulates all coordinates
:param coords: pd.Series of coordinates
:return: boudning list of coordinates
"""
if isinstance(coords, pd.Series):
# assume the series with coords format [(x1, y1), (x2, y2), ...]
coords = coords.to_list()[0]
elif isinstance(coords, pd.DataFrame):
# assume dataframe with X, Y columns
coords = list(zip(coords.X, coords.Y))
xs = [xys[0] for xys in coords]
ys = [xys[1] for xys in coords]
# center
xc = sum(xs) / len(xs)
yc = sum(ys) / len(ys)
# bbox
x1, x2 = min(xs), max(xs)
y1, y2 = min(ys), max(ys)
return (xc, yc), [x1, y1, x2, y2]