From 12070a1acfc5932eee31ef048a78782f698492a4 Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Mon, 17 Dec 2018 14:12:15 +1100 Subject: [PATCH] Plot site .shp as line instead of points --- src/data/csv_to_shp.py | 35 +++++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/src/data/csv_to_shp.py b/src/data/csv_to_shp.py index 46d8f7e..5f8ecc3 100644 --- a/src/data/csv_to_shp.py +++ b/src/data/csv_to_shp.py @@ -1,12 +1,16 @@ """ Converts .csv files to .shape files """ +import os + import click import fiona +import numpy as np import pandas as pd from fiona.crs import from_epsg -from shapely.geometry import Point, mapping +from shapely.geometry import Point, mapping, LineString +from data.parse_shp import convert_coord_systems from utils import setup_logging logger = setup_logging() @@ -25,10 +29,33 @@ def sites_csv_to_shp(input_csv, output_shp): logger.info("Converting %s to %s", input_csv, output_shp) df_sites = pd.read_csv(input_csv, index_col=[0]) logger.info(os.environ.get("GDAL_DATA", None)) - schema = {"geometry": "Point", "properties": {"beach": "str", "site_id": "str"}} + + schema = {"geometry": "LineString", "properties": {"beach": "str", "site_id": "str"}} with fiona.open(output_shp, "w", crs=from_epsg(4326), driver="ESRI Shapefile", schema=schema) as output: for index, row in df_sites.iterrows(): - point = Point(row["x_200_lon"], row["x_200_lat"]) + # Work out where center of profile is + orientation = row["orientation"] + center_profile_x = row["profile_x_lat_lon"] + center_lon = row["lon"] + center_lat = row["lat"] + center_lat_lon = Point(center_lon, center_lat) + center_xy = convert_coord_systems(center_lat_lon) + center_x, center_y = center_xy.xy + + # Work out where landward profile limit is + land_x = center_x + center_profile_x * np.cos(np.deg2rad(orientation)) + land_y = center_y + center_profile_x * np.sin(np.deg2rad(orientation)) + land_xy = Point(land_x, land_y) + land_lat_lon = convert_coord_systems(land_xy, in_coord_system="EPSG:28356", out_coord_system="EPSG:4326") + + # Work out where seaward profile limit is + sea_x = center_x - center_profile_x * np.cos(np.deg2rad(orientation)) + sea_y = center_y - center_profile_x * np.sin(np.deg2rad(orientation)) + sea_xy = Point(sea_x, sea_y) + sea_lat_lon = convert_coord_systems(sea_xy, in_coord_system="EPSG:28356", out_coord_system="EPSG:4326") + + line_string = LineString([land_lat_lon, center_lat_lon, sea_lat_lon]) prop = {"beach": row["beach"], "site_id": index} - output.write({"geometry": mapping(point), "properties": prop}) + output.write({"geometry": mapping(line_string), "properties": prop}) + logger.info("Done!")