Skip to content

Instantly share code, notes, and snippets.

@brey
Created July 13, 2023 09:15
Show Gist options
  • Save brey/667a9372f8a8563deca1cd88e21d3a0b to your computer and use it in GitHub Desktop.
Save brey/667a9372f8a8563deca1cd88e21d3a0b to your computer and use it in GitHub Desktop.
GSHHS
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"id": "47536049",
"metadata": {},
"source": [
"## Prepare latest GSHHS"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69320ecb",
"metadata": {},
"outputs": [],
"source": [
"import geopandas as gp\n",
"import numpy as np\n",
"import pandas as pd\n",
"import shapely"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ffcdc723",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib as mpl\n",
"mpl.rcParams['agg.path.chunksize'] = 10000"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69687495",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aff3260c",
"metadata": {},
"outputs": [],
"source": [
"import cartopy.feature as cf"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7b3d2407",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget"
]
},
{
"cell_type": "markdown",
"id": "73313b40",
"metadata": {},
"source": [
"#### Old way"
]
},
{
"cell_type": "markdown",
"id": "2d176eb0",
"metadata": {},
"source": [
"Using the cartopy was working but then the format changed "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9c2fe320",
"metadata": {},
"outputs": [],
"source": [
"gi = cf.GSHHSFeature(scale=\"full\", levels=[1])\n",
"wo = gp.GeoDataFrame(geometry=[x for x in gi.geometries()])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c42f73bc",
"metadata": {},
"outputs": [],
"source": [
"wo = wo.explode(index_parts=True).reset_index(drop=True)\n",
"wo['area'] = wo.area\n",
"wo = wo.sort_values(by=\"area\", ascending=0) # sort\n",
"wo = wo.reset_index(drop=True)\n",
"#wo_ = wo.iloc[:100]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "928517ef",
"metadata": {},
"outputs": [],
"source": [
"wo_.plot()"
]
},
{
"cell_type": "markdown",
"id": "bd391394",
"metadata": {},
"source": [
"#### New way"
]
},
{
"cell_type": "markdown",
"id": "2c30e64d",
"metadata": {},
"source": [
"Download the latest shapefile : http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d71fa32e",
"metadata": {},
"outputs": [],
"source": [
"w1 = gp.read_file('/Users/brey/.local/share/cartopy/shapefiles/gshhg-shp-2/GSHHS_shp/f/GSHHS_f_L1.shp')\n",
"w2 = gp.read_file('/Users/brey/.local/share/cartopy/shapefiles/gshhg-shp-2/GSHHS_shp/f/GSHHS_f_L6.shp')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "45785acb",
"metadata": {},
"outputs": [],
"source": [
"w5 = gp.read_file('/Users/brey/.local/share/cartopy/shapefiles/gshhg-shp-2/GSHHS_shp/f/GSHHS_f_L5.shp')"
]
},
{
"cell_type": "markdown",
"id": "fc76412a",
"metadata": {},
"source": [
"First check fails"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cb564170",
"metadata": {},
"outputs": [],
"source": [
"w5.buffer(0).is_valid.all(), (w5.buffer(0).boundary.geom_type == \"LineString\").all()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "67dbfcc1",
"metadata": {},
"outputs": [],
"source": [
"# fix the first polygon\n",
"x, y = w5.iloc[[0]].boundary[0].coords.xy\n",
"\n",
"x = x[2:-1]\n",
"y = y[2:-1]\n",
"\n",
"\n",
"x.append(x[-1])\n",
"y.append(-90.)\n",
"x.append(x[0])\n",
"y.append(y[0])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8ea33bde",
"metadata": {},
"outputs": [],
"source": [
"w5.loc[0:0, \"geometry\"] = shapely.Polygon(zip(x,y))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6917a4b2",
"metadata": {},
"outputs": [],
"source": [
"x, y = w5.iloc[[1]].boundary[1].coords.xy\n",
"x[0] = 0 \n",
"x[1] = 0\n",
"x[-2] = 0 \n",
"x[-1] = 0"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "feacc558",
"metadata": {},
"outputs": [],
"source": [
"w5.loc[1:1, \"geometry\"] = shapely.Polygon(zip(x,y))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "72046e54",
"metadata": {},
"outputs": [],
"source": [
"w5_ = gp.GeoDataFrame(geometry=w5.buffer(0))"
]
},
{
"cell_type": "markdown",
"id": "bb0f7528",
"metadata": {},
"source": [
"After the fix, it still fails"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3bfdcc8f",
"metadata": {},
"outputs": [],
"source": [
"w5_.buffer(0).is_valid.all(), (w5_.buffer(0).boundary.geom_type == \"LineString\").all()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d6c2b61d",
"metadata": {},
"outputs": [],
"source": [
"w5_.iloc[[0]].boundary.explode(index_parts=True).plot(cmap='OrRd')\n",
"s1=plt.gca()"
]
},
{
"cell_type": "markdown",
"id": "910db17a",
"metadata": {},
"source": [
"### merge"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d988c5b3",
"metadata": {},
"outputs": [],
"source": [
"w1_ = gp.GeoDataFrame(geometry=w1.buffer(0))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b0f5edd1",
"metadata": {},
"outputs": [],
"source": [
"w = gp.pd.concat([w1_,w5_])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "05599869",
"metadata": {},
"outputs": [],
"source": [
"cols = [a for a in w.columns if a not in [\"geometry\",\"area\"]]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "47da32cb",
"metadata": {},
"outputs": [],
"source": [
"w = w.drop(columns=cols)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe99fb36",
"metadata": {},
"outputs": [],
"source": [
"w['area'] = w.to_crs(epsg=4087).area\n",
"w = w.explode(index_parts=True).reset_index(drop=True)\n",
"w = w.sort_values(by=\"area\", ascending=0) # sort\n",
"w = w.reset_index(drop=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "000ad520",
"metadata": {},
"outputs": [],
"source": [
"w.boundary.plot()"
]
},
{
"cell_type": "markdown",
"id": "15348cd6",
"metadata": {},
"source": [
"#### Simplify"
]
},
{
"cell_type": "markdown",
"id": "a6a02703",
"metadata": {},
"source": [
"The `pyposeidon` utility can handle that input"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7401b447",
"metadata": {},
"outputs": [],
"source": [
"from pyposeidon.utils.coastfix import simplify"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3ac3831e",
"metadata": {},
"outputs": [],
"source": [
"coasts = simplify(w)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b7d4f26f",
"metadata": {},
"outputs": [],
"source": [
"coasts.is_valid.all()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e5d1dc03",
"metadata": {},
"outputs": [],
"source": [
"(coasts.boundary.geom_type == \"LineString\").all() # no multiLines"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a23456d4",
"metadata": {},
"outputs": [],
"source": [
"coasts.boundary.plot()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b18628e2",
"metadata": {},
"outputs": [],
"source": [
"#### save to file\n",
"coasts.to_file('./coasts')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "dev",
"language": "python",
"name": "dev"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment