Skip to content

Instantly share code, notes, and snippets.

@OlegJakushkin

OlegJakushkin/Rome.ipynb

Last active February 12, 2021 11:19
Show Gist options
  • Save OlegJakushkin/506277ed98dd6a7a79c9082a4a158ecd to your computer and use it in GitHub Desktop.
Save OlegJakushkin/506277ed98dd6a7a79c9082a4a158ecd to your computer and use it in GitHub Desktop.
Download ZIP
olejak/gis
Raw
geotiff-to-geojson.sh
#!/bin/bash
#
# Requires:
# - gdal_sieve.py
# - ogr2ogr (GDAL)
# - topojson (node.js)
# Grab the relative directory for source file.
SRC_DIR=`dirname $0`
# Which raster to compress.
ORG_FILE="$SRC_DIR/Rome-DEM.tif"
# Final output file.
OUTPUT_FILE="$SRC_DIR/Rome.json"
echo "Processing $ORG_FILE."
# Where to output the new file.
TMP_DIR=./tmp
# The amount of times the file should be passed over.
ITERATIONS=3
# Threshold for each iteration.
THRESHOLD=40
# TopoJSON area threshold for simplification.
TOPO_COMPRESSION=0.005
# Setup internal vars.
_CUR=$THRESHOLD
_COMPRESSION=$(($ITERATIONS * $THRESHOLD))
rm -rf $TMP_DIR
mkdir -p $TMP_DIR
# Start sieve passes.
gdal_sieve.py -st $THRESHOLD -4 $ORG_FILE $TMP_DIR/output-"$THRESHOLD".tiff
while [ $_CUR -le $_COMPRESSION ]; do
let _PREV=$_CUR
let _CUR=$_CUR+$THRESHOLD
echo "Compressing output-$_PREV.tiff into $_CUR.tiff"
gdal_sieve.py -st $THRESHOLD -4 "$TMP_DIR/output-$_PREV.tiff" \
"$TMP_DIR/output-$_CUR.tiff"
rm "$TMP_DIR/output-$_PREV.tiff"
done
# Raster to vector.
gdal_polygonize.py $TMP_DIR/output-"$_CUR".tiff \
-f "ESRI Shapefile" $TMP_DIR vector n
# Change shapefile to geojson without the 0 layer, which is water.
ogr2ogr -f "GeoJSON" -where "n != 0" "$TMP_DIR/geojson.json" $TMP_DIR/vector.shp
# Convert to compressed TopoJSON.
geo2topo -o $OUTPUT_FILE \
--no-stitch-poles \
-s $TOPO_COMPRESSION \
-p -- "$TMP_DIR/geojson.json"
# Clean up.
rm -rf $TMP_DIR
Display the source blob
Display the rendered blob
Raw
Rome.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "current-cannon",
"metadata": {},
"outputs": [],
"source": [
"# docker run --rm -p 10000:8888 -e JUPYTER_ENABLE_LAB=yes --user root -e NB_GID=100 -e GEN_CERT=yes -e GRANT_SUDO=yes -v C:\\TestDir:/home/jovyan/work --name ojgis2 jupyter/datascience-notebook:lab-2.2.9\n",
"from jupyter/datascience-notebook:lab-2.2.9\n",
"USER 0\n",
"ARG MapboxAccessTokenJupyter=pk.eyJ1Ijoic3VwZXJpb3IwIiwiYSI6ImNrbDF3NThvaTAyZGgyb3J3OGpiaGx0Y3IifQ.5_I3PuByqUc-PkDBMuSzUQ\n",
"\n",
"RUN apt-get update && \\\n",
" apt-get install -y software-properties-common python-numpy htop mc gdal-bin libgdal-dev && \\\n",
" pip install keplergl elevation rasterio topojson && \\\n",
" pip uninstall -y keplergl\n",
"\n",
"RUN npm install -g --force yarn topojson-server topojson && \\\n",
" pip install GDAL==3.0.4 && \\\n",
" pip install raster2xyz\n",
"\n",
"RUN echo \"fs.inotify.max_user_watches=524288\" >> /etc/sysctl.conf\n",
"\n",
"RUN git clone --recursive https://github.com/keplergl/kepler.gl && \\\n",
" cd ./kepler.gl/bindings/kepler.gl-jupyter && cd js && yarn && npm run build && \\\n",
" cd ../ && pip install -e .\n",
"\n",
"RUN cd ./kepler.gl/bindings/kepler.gl-jupyter \\\n",
" && jupyter nbextension install --py --symlink --sys-prefix keplergl \\\n",
" && jupyter nbextension enable --py --sys-prefix keplergl\n",
"\n",
"RUN cd ./kepler.gl/bindings/kepler.gl-jupyter \\\n",
" && jupyter labextension install @jupyter-widgets/jupyterlab-manager \\\n",
" && jupyter labextension install js"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "solved-scoop",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "thermal-hanging",
"metadata": {},
"outputs": [],
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "sexual-birth",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"'GDAL_DATA' in os.environ\n",
"!fio env --gdal-data\n",
"!export GDAL_DATA=/opt/conda/lib/python3.8/site-packages/fiona/gdal_data\n",
"os.environ['GDAL_DATA'] = '/opt/conda/lib/python3.8/site-packages/fiona/gdal_data'\n",
"os.environ['MapboxAccessTokenJupyter'] = 'pk.eyJ1Ijoic3VwZXJpb3IwIiwiYSI6ImNrbDF3NThvaTAyZGgyb3J3OGpiaGx0Y3IifQ.5_I3PuByqUc-PkDBMuSzUQ'"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "czech-relevance",
"metadata": {},
"outputs": [],
"source": [
"!eio selfcheck\n",
"!eio clip -o Rome-DEM.tif --bounds 12.35 41.8 12.65 42"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "falling-graham",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "coordinate-moment",
"metadata": {},
"outputs": [],
"source": [
"import rasterio\n",
"import rasterio.features\n",
"import rasterio.warp\n",
"from rasterio.plot import show\n",
"from matplotlib import pyplot\n",
"with rasterio.open('Rome-DEM.tif') as dataset:\n",
"\n",
" # Read the dataset's valid data mask as a ndarray.\n",
" mask = dataset.dataset_mask()\n",
"\n",
" # Extract feature shapes and values from the array.\n",
" for geom, val in rasterio.features.shapes(\n",
" mask, transform=dataset.transform):\n",
"\n",
" # Transform shapes from the dataset's own coordinate\n",
" # reference system to CRS84 (EPSG:4326).\n",
" geom = rasterio.warp.transform_geom(\n",
" dataset.crs, 'EPSG:4326', geom, precision=6)\n",
" \n",
" fig, ax = pyplot.subplots(1, figsize=(12, 12))\n",
" #show(dataset.read(), transform=dataset.transform)\n",
" show((dataset, 1), cmap='Greys_r', interpolation='none', ax=ax)\n",
" show((dataset, 1), contour=True, ax=ax)\n",
" pyplot.show()\n",
"!rio shapes 'Rome-DEM.tif' --bidx 1 --precision 6 > shade.geojson"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "helpful-incident",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "coordinate-marketplace",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "twelve-despite",
"metadata": {},
"outputs": [],
"source": [
"!./geotiff-to-geojson.sh"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "mysterious-feedback",
"metadata": {},
"outputs": [],
"source": [
"from raster2xyz.raster2xyz import Raster2xyz\n",
"\n",
"input_raster = \"Rome-DEM.tif\"\n",
"out_csv = \"Rome.csv\"\n",
"\n",
"rtxyz = Raster2xyz()\n",
"rtxyz.translate(input_raster, out_csv)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "concrete-fancy",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"User Guide: https://docs.kepler.gl/docs/keplergl-jupyter\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "761cd593fd674ea5a643867222004055",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"KeplerGl(config={'version': 'v1', 'config': {'visState': {'filters': [], 'layers': [{'id': 'xxu0rsf', 'type': …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from keplergl import KeplerGl\n",
"import pandas as pd\n",
"from geopandas import GeoDataFrame\n",
"from shapely.geometry import Point\n",
"import fiona\n",
"\n",
"\n",
"df = pd.read_csv('Rome.csv')\n",
"\n",
"geometry = [Point(xy) for xy in zip(df.x, df.y)]\n",
"crs = {'init': 'epsg:2263'} #http://www.spatialreference.org/ref/epsg/2263/\n",
"geo_df = GeoDataFrame(df, crs=crs, geometry=geometry)\n",
"config = {\n",
" \"version\": \"v1\",\n",
" \"config\": {\n",
" \"visState\": {\n",
" \"filters\": [],\n",
" \"layers\": [\n",
" {\n",
" \"id\": \"xxu0rsf\",\n",
" \"type\": \"point\",\n",
" \"config\": {\n",
" \"dataId\": \"depth\",\n",
" \"label\": \"depth\",\n",
" \"color\": [\n",
" 18,\n",
" 147,\n",
" 154\n",
" ],\n",
" \"columns\": {\n",
" \"lat\": \"y\",\n",
" \"lng\": \"x\",\n",
" \"altitude\": \"z\"\n",
" },\n",
" \"isVisible\": True,\n",
" \"visConfig\": {\n",
" \"radius\": 5.2,\n",
" \"fixedRadius\": False,\n",
" \"opacity\": 0.59,\n",
" \"outline\": False,\n",
" \"thickness\": 0.5,\n",
" \"strokeColor\": None,\n",
" \"colorRange\": {\n",
" \"name\": \"Global Warming\",\n",
" \"type\": \"sequential\",\n",
" \"category\": \"Uber\",\n",
" \"colors\": [\n",
" \"#5A1846\",\n",
" \"#900C3F\",\n",
" \"#C70039\",\n",
" \"#E3611C\",\n",
" \"#F1920E\",\n",
" \"#FFC300\"\n",
" ]\n",
" },\n",
" \"strokeColorRange\": {\n",
" \"name\": \"Global Warming\",\n",
" \"type\": \"sequential\",\n",
" \"category\": \"Uber\",\n",
" \"colors\": [\n",
" \"#5A1846\",\n",
" \"#900C3F\",\n",
" \"#C70039\",\n",
" \"#E3611C\",\n",
" \"#F1920E\",\n",
" \"#FFC300\"\n",
" ]\n",
" },\n",
" \"radiusRange\": [\n",
" 0,\n",
" 50\n",
" ],\n",
" \"filled\": True\n",
" },\n",
" \"hidden\": False,\n",
" \"textLabel\": [\n",
" {\n",
" \"field\": None,\n",
" \"color\": [\n",
" 255,\n",
" 255,\n",
" 255\n",
" ],\n",
" \"size\": 18,\n",
" \"offset\": [\n",
" 0,\n",
" 0\n",
" ],\n",
" \"anchor\": \"start\",\n",
" \"alignment\": \"center\"\n",
" }\n",
" ]\n",
" },\n",
" \"visualChannels\": {\n",
" \"colorField\": {\n",
" \"name\": \"z\",\n",
" \"type\": \"integer\"\n",
" },\n",
" \"colorScale\": \"quantile\",\n",
" \"strokeColorField\": None,\n",
" \"strokeColorScale\": \"quantile\",\n",
" \"sizeField\": None,\n",
" \"sizeScale\": \"linear\"\n",
" }\n",
" }\n",
" ],\n",
" \"interactionConfig\": {\n",
" \"tooltip\": {\n",
" \"fieldsToShow\": {\n",
" \"depth\": [\n",
" {\n",
" \"name\": \"x\",\n",
" \"format\": None\n",
" },\n",
" {\n",
" \"name\": \"y\",\n",
" \"format\": None\n",
" },\n",
" {\n",
" \"name\": \"z\",\n",
" \"format\": None\n",
" }\n",
" ]\n",
" },\n",
" \"compareMode\": False,\n",
" \"compareType\": \"absolute\",\n",
" \"enabled\": True\n",
" },\n",
" \"brush\": {\n",
" \"size\": 0.5,\n",
" \"enabled\": False\n",
" },\n",
" \"geocoder\": {\n",
" \"enabled\": False\n",
" },\n",
" \"coordinate\": {\n",
" \"enabled\": False\n",
" }\n",
" },\n",
" \"layerBlending\": \"normal\",\n",
" \"splitMaps\": [],\n",
" \"animationConfig\": {\n",
" \"currentTime\": None,\n",
" \"speed\": 1\n",
" }\n",
" },\n",
" \"mapState\": {\n",
" \"bearing\": 12.5817852288174,\n",
" \"dragRotate\": True,\n",
" \"latitude\": 41.89028347702737,\n",
" \"longitude\": 12.462709954174432,\n",
" \"pitch\": 59.50249662967758,\n",
" \"zoom\": 13.769164917078767,\n",
" \"isSplit\": False\n",
" },\n",
" \"mapStyle\": {\n",
" \"styleType\": \"dark\",\n",
" \"topLayerGroups\": {},\n",
" \"visibleLayerGroups\": {\n",
" \"label\": True,\n",
" \"road\": True,\n",
" \"border\": False,\n",
" \"building\": True,\n",
" \"water\": True,\n",
" \"land\": True,\n",
" \"3d building\": False\n",
" },\n",
" \"threeDBuildingColor\": [\n",
" 9.665468314072013,\n",
" 17.18305478057247,\n",
" 31.1442867897876\n",
" ],\n",
" \"mapStyles\": {}\n",
" }\n",
" }\n",
"}\n",
"\n",
"map_1 = KeplerGl(height=800, config=config)\n",
"map_1.add_data(data=geo_df, name='depth')\n",
"map_1"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.8.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
@OlegJakushkin
Copy link
Author

Shows:
image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment