{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### `geom_imshow()` renders `NaN` values as transparent pixels"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import rasterio as ro\n",
"\n",
"from lets_plot import *\n",
"LetsPlot.setup_html()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Load Data - Japan Digital Elevation Model (DEM)\n",
"\n",
"Data are available at\n",
"\n",
"[Harvard Dataverse](https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/28762)\n",
"\n",
"Japan DEM (c) USGS. Processed for Japan by Skinner Regional Systems Analysis. Cambridge: Center for Geographic Analysis, Harvard University, 2015."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Load data with Rasterio\n",
"\n",
"jp_tif=ro.open('./data/japan_dem_wgs84/japan_dem_wgs84.tif')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EPSG:4326\n",
"BoundingBox(left=127.9421598662918, bottom=29.201489591057303, right=155.88170436191123, top=50.02079185141284)\n"
]
}
],
"source": [
"# Metadata\n",
"\n",
"jp_crs=jp_tif.crs\n",
"jp_bounds=jp_tif.bounds\n",
"\n",
"print(\"{}\\n{}\".format(jp_crs, jp_bounds))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# create the \"extent\" to use in geom_imshow()\n",
"\n",
"xmin=jp_bounds[0]\n",
"xmax=jp_bounds[2]\n",
"ymin=jp_bounds[1]\n",
"ymax=jp_bounds[3]\n",
"\n",
"jp_ext=[xmin,xmax, ymin, ymax]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Get the first band as a 2D numpy array\n",
"\n",
"jp_arr=jp_tif.read(1)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(numpy.ndarray, (2038, 2735))"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(jp_arr), np.shape(jp_arr)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ggplot() + geom_imshow(jp_arr, extent=jp_ext)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Easy to see (chart below) that values in `jp_arr` form 2 clusters: large negative values and values close to 0. \n",
"\n",
"After \"normalization\", those 2 clasters fall into opposite sides of the luminosity spectrum, and as result we can only see blak&wite picture."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(ggplot() \n",
" + geom_histogram(aes(x=jp_arr.flatten()))\n",
" + geom_vline(xintercept=-32768, color=\"red\")\n",
" + geom_label(label=\"-32768\", x=-30000, y=5.2E6)\n",
" + ggtitle(\"Distribution of elevations in Japan DEM\")\n",
" + xlab(\"elevation (m)\")\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Re-encode the sea level using `NaN` values\n",
"\n",
"The sea level elevation is traditionally encoded as minimal 16-bit value: -32768. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"jp_arr_nan=jp_arr.copy().astype(np.float32)\n",
"jp_arr_nan[jp_arr<-30000]=np.nan\n",
"\n",
"(ggplot() \n",
" + geom_histogram(aes(x=jp_arr_nan.flatten()))\n",
" + ggtitle(\"Distribution of elevations in Japan DEM\", subtitle=\"excluding 'sea level'\")\n",
" + xlab(\"elevation (m)\")\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sea level is now transparent"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ggplot() + geom_imshow(jp_arr_nan, extent=jp_ext)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'xmin': 127.9421598662918,\n",
" 'ymin': 29.201489591057303,\n",
" 'xmax': 155.88170436191123,\n",
" 'ymax': 50.02079185141284}"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"jp_bbox = dict(zip(['xmin', 'ymin', 'xmax', 'ymax'], jp_bounds))\n",
"jp_bbox"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"labels = dict(\n",
" label=[\"Sea of Japan\", \"Pacific\\nOcean\"],\n",
" x=[134, 144],\n",
" y=[41, 33]\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(ggplot() \n",
" + geom_rect(**jp_bbox, fill='#aadaff', alpha=0.2)\n",
" + geom_text(aes(\"x\", \"y\", label=\"label\"), data=labels, color='#578bcc', size=11, fontface='italic', family=\"Courier\")\n",
" + geom_imshow(jp_arr_nan, \"viridis\", extent=jp_ext, vmax=1500)\n",
" + coord_fixed(xlim=[128,147], ylim=[30, 46])\n",
" + theme_bw() \n",
" + theme(axis_title=\"blank\", plot_title=element_text(face=\"bold\"))\n",
" + ggsize(800, 680)\n",
" + ggtitle(\"Japan Digital Elevation Model\") \n",
" + labs(caption=\"Japan DEM (c) USGS.\\n\" +\n",
" \"Processed for Japan by Skinner Regional Systems Analysis.\\n\" +\n",
" \"Cambridge: Center for Geographic Analysis, Harvard University, 2015.\")\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.15"
}
},
"nbformat": 4,
"nbformat_minor": 4
}