{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### `geom_imshow()` renders `NaN` values as transparent pixels"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:24.202251Z",
"iopub.status.busy": "2025-11-05T13:50:24.202142Z",
"iopub.status.idle": "2025-11-05T13:50:24.227639Z",
"shell.execute_reply": "2025-11-05T13:50:24.227341Z"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"from osgeo import gdal, osr\n",
"\n",
"from lets_plot import *"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:24.229005Z",
"iopub.status.busy": "2025-11-05T13:50:24.228922Z",
"iopub.status.idle": "2025-11-05T13:50:24.230712Z",
"shell.execute_reply": "2025-11-05T13:50:24.230551Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"LetsPlot.setup_html()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:24.246046Z",
"iopub.status.busy": "2025-11-05T13:50:24.245953Z",
"iopub.status.idle": "2025-11-05T13:50:24.248117Z",
"shell.execute_reply": "2025-11-05T13:50:24.247757Z"
}
},
"outputs": [],
"source": [
"def get_bounds(data):\n",
" xmin, xpixel, _, ymax, _, ypixel = data.GetGeoTransform()\n",
" width, height = data.RasterXSize, data.RasterYSize\n",
" xmax = xmin + width * xpixel\n",
" ymin = ymax + height * ypixel\n",
" return xmin, ymin, xmax, ymax\n",
"\n",
"def get_crs(data):\n",
" proj = osr.SpatialReference(wkt=data.GetProjection())\n",
" return \"EPSG:{0}\".format(proj.GetAttrValue('AUTHORITY', 1))"
]
},
{
"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": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:24.249139Z",
"iopub.status.busy": "2025-11-05T13:50:24.249055Z",
"iopub.status.idle": "2025-11-05T13:50:25.513300Z",
"shell.execute_reply": "2025-11-05T13:50:25.512795Z"
}
},
"outputs": [],
"source": [
"# Load data with GDAL\n",
"\n",
"jp_tif = gdal.Open(\"https://raw.githubusercontent.com/JetBrains/lets-plot-docs/master/data/japan_dem_wgs84/japan_dem_wgs84.tif\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:25.514718Z",
"iopub.status.busy": "2025-11-05T13:50:25.514592Z",
"iopub.status.idle": "2025-11-05T13:50:25.531735Z",
"shell.execute_reply": "2025-11-05T13:50:25.531477Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EPSG:4326\n",
"(127.9421598662918, 29.201489591057303, 155.88170436191123, 50.02079185141284)\n"
]
}
],
"source": [
"# Metadata\n",
"\n",
"jp_crs = get_crs(jp_tif)\n",
"jp_bounds = get_bounds(jp_tif)\n",
"\n",
"print(\"{}\\n{}\".format(jp_crs, jp_bounds))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:25.532669Z",
"iopub.status.busy": "2025-11-05T13:50:25.532597Z",
"iopub.status.idle": "2025-11-05T13:50:25.534014Z",
"shell.execute_reply": "2025-11-05T13:50:25.533835Z"
}
},
"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": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:25.534913Z",
"iopub.status.busy": "2025-11-05T13:50:25.534842Z",
"iopub.status.idle": "2025-11-05T13:50:25.542157Z",
"shell.execute_reply": "2025-11-05T13:50:25.541952Z"
}
},
"outputs": [],
"source": [
"# Get the first band as a 2D numpy array\n",
"\n",
"jp_arr = jp_tif.GetRasterBand(1).ReadAsArray()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:25.543172Z",
"iopub.status.busy": "2025-11-05T13:50:25.543085Z",
"iopub.status.idle": "2025-11-05T13:50:25.545245Z",
"shell.execute_reply": "2025-11-05T13:50:25.545055Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(numpy.ndarray, (2038, 2735))"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(jp_arr), np.shape(jp_arr)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:25.546076Z",
"iopub.status.busy": "2025-11-05T13:50:25.545997Z",
"iopub.status.idle": "2025-11-05T13:50:25.730470Z",
"shell.execute_reply": "2025-11-05T13:50:25.730047Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 9,
"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": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:25.732040Z",
"iopub.status.busy": "2025-11-05T13:50:25.731939Z",
"iopub.status.idle": "2025-11-05T13:50:30.954237Z",
"shell.execute_reply": "2025-11-05T13:50:30.953879Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 10,
"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": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:30.955414Z",
"iopub.status.busy": "2025-11-05T13:50:30.955306Z",
"iopub.status.idle": "2025-11-05T13:50:33.795731Z",
"shell.execute_reply": "2025-11-05T13:50:33.795323Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 11,
"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": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:33.796693Z",
"iopub.status.busy": "2025-11-05T13:50:33.796594Z",
"iopub.status.idle": "2025-11-05T13:50:33.998093Z",
"shell.execute_reply": "2025-11-05T13:50:33.997784Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ggplot() + geom_imshow(jp_arr_nan, extent=jp_ext)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:34.000904Z",
"iopub.status.busy": "2025-11-05T13:50:34.000827Z",
"iopub.status.idle": "2025-11-05T13:50:34.002875Z",
"shell.execute_reply": "2025-11-05T13:50:34.002655Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'xmin': 127.9421598662918,\n",
" 'ymin': 29.201489591057303,\n",
" 'xmax': 155.88170436191123,\n",
" 'ymax': 50.02079185141284}"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"jp_bbox = dict(zip(['xmin', 'ymin', 'xmax', 'ymax'], jp_bounds))\n",
"jp_bbox"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:34.003701Z",
"iopub.status.busy": "2025-11-05T13:50:34.003626Z",
"iopub.status.idle": "2025-11-05T13:50:34.004941Z",
"shell.execute_reply": "2025-11-05T13:50:34.004775Z"
}
},
"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": 15,
"metadata": {
"execution": {
"iopub.execute_input": "2025-11-05T13:50:34.005730Z",
"iopub.status.busy": "2025-11-05T13:50:34.005659Z",
"iopub.status.idle": "2025-11-05T13:50:34.189403Z",
"shell.execute_reply": "2025-11-05T13:50:34.189140Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 15,
"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.10.13"
}
},
"nbformat": 4,
"nbformat_minor": 4
}