{
"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": "2024-04-17T07:36:22.794779Z",
"iopub.status.busy": "2024-04-17T07:36:22.794695Z",
"iopub.status.idle": "2024-04-17T07:36:23.130487Z",
"shell.execute_reply": "2024-04-17T07:36:23.130144Z"
}
},
"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": "2024-04-17T07:36:23.131885Z",
"iopub.status.busy": "2024-04-17T07:36:23.131754Z",
"iopub.status.idle": "2024-04-17T07:36:23.133860Z",
"shell.execute_reply": "2024-04-17T07:36:23.133689Z"
}
},
"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": "2024-04-17T07:36:23.146965Z",
"iopub.status.busy": "2024-04-17T07:36:23.146798Z",
"iopub.status.idle": "2024-04-17T07:36:23.148721Z",
"shell.execute_reply": "2024-04-17T07:36:23.148519Z"
}
},
"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": "2024-04-17T07:36:23.149641Z",
"iopub.status.busy": "2024-04-17T07:36:23.149569Z",
"iopub.status.idle": "2024-04-17T07:36:24.446485Z",
"shell.execute_reply": "2024-04-17T07:36:24.446126Z"
}
},
"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": "2024-04-17T07:36:24.447848Z",
"iopub.status.busy": "2024-04-17T07:36:24.447751Z",
"iopub.status.idle": "2024-04-17T07:36:24.464521Z",
"shell.execute_reply": "2024-04-17T07:36:24.464222Z"
}
},
"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": "2024-04-17T07:36:24.465515Z",
"iopub.status.busy": "2024-04-17T07:36:24.465424Z",
"iopub.status.idle": "2024-04-17T07:36:24.467270Z",
"shell.execute_reply": "2024-04-17T07:36:24.467092Z"
}
},
"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": "2024-04-17T07:36:24.468246Z",
"iopub.status.busy": "2024-04-17T07:36:24.468103Z",
"iopub.status.idle": "2024-04-17T07:36:24.474589Z",
"shell.execute_reply": "2024-04-17T07:36:24.474378Z"
}
},
"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": "2024-04-17T07:36:24.475705Z",
"iopub.status.busy": "2024-04-17T07:36:24.475548Z",
"iopub.status.idle": "2024-04-17T07:36:24.478120Z",
"shell.execute_reply": "2024-04-17T07:36:24.477948Z"
}
},
"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": "2024-04-17T07:36:24.479105Z",
"iopub.status.busy": "2024-04-17T07:36:24.478991Z",
"iopub.status.idle": "2024-04-17T07:36:24.560702Z",
"shell.execute_reply": "2024-04-17T07:36:24.560458Z"
}
},
"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": "2024-04-17T07:36:24.562123Z",
"iopub.status.busy": "2024-04-17T07:36:24.562018Z",
"iopub.status.idle": "2024-04-17T07:36:28.617525Z",
"shell.execute_reply": "2024-04-17T07:36:28.617168Z"
}
},
"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": "2024-04-17T07:36:28.619311Z",
"iopub.status.busy": "2024-04-17T07:36:28.619196Z",
"iopub.status.idle": "2024-04-17T07:36:31.121906Z",
"shell.execute_reply": "2024-04-17T07:36:31.121462Z"
}
},
"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": "2024-04-17T07:36:31.123750Z",
"iopub.status.busy": "2024-04-17T07:36:31.123655Z",
"iopub.status.idle": "2024-04-17T07:36:31.236832Z",
"shell.execute_reply": "2024-04-17T07:36:31.236360Z"
}
},
"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": "2024-04-17T07:36:31.240370Z",
"iopub.status.busy": "2024-04-17T07:36:31.240250Z",
"iopub.status.idle": "2024-04-17T07:36:31.242674Z",
"shell.execute_reply": "2024-04-17T07:36:31.242498Z"
}
},
"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": "2024-04-17T07:36:31.244146Z",
"iopub.status.busy": "2024-04-17T07:36:31.243980Z",
"iopub.status.idle": "2024-04-17T07:36:31.245458Z",
"shell.execute_reply": "2024-04-17T07:36:31.245282Z"
}
},
"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": "2024-04-17T07:36:31.246952Z",
"iopub.status.busy": "2024-04-17T07:36:31.246840Z",
"iopub.status.idle": "2024-04-17T07:36:31.338132Z",
"shell.execute_reply": "2024-04-17T07:36:31.337700Z"
}
},
"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
}