reproject knmi radar using rasterio.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "given following packages\n",
    "---\n",
    "- `gdal` `2.4.4`\n",
    "- `rasterio` `1.1.5`\n",
    "- `proj 7.0.0`\n",
    "- `pyproj` `2.6.1.post1`\n",
    "- `h5py` `2.10.0`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import h5py\n",
    "import rasterio\n",
    "from rasterio.transform import from_origin\n",
    "from rasterio.warp import calculate_default_transform, reproject, Resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load KNMI radar data\n",
    "filename_h5 = r'RAD_NL25_PCP_NA_201910281900.h5'\n",
    "\n",
    "# Filename h5 converted to tif\n",
    "filename_rad = 'RAD_NL25_PCP_NA_201910281900_rad.tif'\n",
    "\n",
    "# Filename tif reprojected to RD\n",
    "filename_rd = 'RAD_NL25_PCP_NA_201910281900_rd_rio.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open h5 file into numpy array\n",
    "f = h5py.File(filename_h5, 'r')\n",
    "im = np.copy(f['image1']['image_data'])\n",
    "im = im.astype(np.float32)\n",
    "f.close()\n",
    "arr = np.ma.masked_where(im==255., im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save array to tif\n",
    "transform = from_origin(0, -3650, 1, 1)\n",
    "new_dataset = rasterio.open(filename_rad, 'w', driver='GTiff',\n",
    "                            height = arr.shape[0], width = arr.shape[1],\n",
    "                            count=1, dtype=str(arr.dtype),\n",
    "                            crs='+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752',                            \n",
    "                            transform=transform)\n",
    "\n",
    "new_dataset.write(arr, 1)\n",
    "new_dataset.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reproject to RD projection\n",
    "dst_crs = rasterio.crs.CRS.from_proj4('+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +<>')\n",
    "#dst_crs = rasterio.crs.CRS.from_epsg(4326)\n",
    "\n",
    "with rasterio.open(filename_rad) as src:\n",
    "    transform, width, height = calculate_default_transform(\n",
    "        src.crs, dst_crs, src.width, src.height, *src.bounds)\n",
    "    kwargs = src.meta.copy()\n",
    "    kwargs.update({\n",
    "        'crs': dst_crs,\n",
    "        'transform': transform,\n",
    "        'width': width,\n",
    "        'height': height\n",
    "    })\n",
    "\n",
    "    with rasterio.open(filename_rd, 'w', **kwargs) as dst:\n",
    "        for i in range(1, src.count + 1):\n",
    "            reproject(\n",
    "                source=rasterio.band(src, i),\n",
    "                destination=rasterio.band(dst, i),\n",
    "                src_transform=src.transform,\n",
    "                src_crs=src.crs,\n",
    "                dst_transform=transform,\n",
    "                dst_crs=dst_crs,\n",
    "                resampling=Resampling.nearest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Driver: GTiff/GeoTIFF\n",
      "Files: RAD_NL25_PCP_NA_201910281900_rd_rio.tif\n",
      "Size is 768, 826\n",
      "Coordinate System is:\n",
      "PROJCS[\"unnamed\",\n",
      "    GEOGCS[\"Bessel 1841\",\n",
      "        DATUM[\"unknown\",\n",
      "            SPHEROID[\"bessel\",6377397.155,299.1528128]],\n",
      "        PRIMEM[\"Greenwich\",0],\n",
      "        UNIT[\"degree\",0.0174532925199433]],\n",
      "    PROJECTION[\"Oblique_Stereographic\"],\n",
      "    PARAMETER[\"latitude_of_origin\",52.15616055555555],\n",
      "    PARAMETER[\"central_meridian\",5.38763888888889],\n",
      "    PARAMETER[\"scale_factor\",0.9999079],\n",
      "    PARAMETER[\"false_easting\",155000],\n",
      "    PARAMETER[\"false_northing\",463000],\n",
      "    UNIT[\"metre\",1,\n",
      "        AUTHORITY[\"EPSG\",\"9001\"]]]\n",
      "Origin = (-236275.402756847790442,900792.860893698176369)\n",
      "Pixel Size = (961.025429589943087,-961.025429589943087)\n",
      "Metadata:\n",
      "  AREA_OR_POINT=Area\n",
      "Image Structure Metadata:\n",
      "  INTERLEAVE=BAND\n",
      "Corner Coordinates:\n",
      "Upper Left  ( -236275.403,  900792.861) (  0d52'23.92\"W, 55d55'59.21\"N)\n",
      "Lower Left  ( -236275.403,  106985.856) (  0d 3'31.53\"E, 48d49'46.56\"N)\n",
      "Upper Right (  501792.127,  900792.861) ( 10d56'25.11\"E, 55d57'59.00\"N)\n",
      "Lower Right (  501792.127,  106985.856) ( 10d 6'45.81\"E, 48d51'23.91\"N)\n",
      "Center      (  132758.362,  503889.358) (  5d 3'35.58\"E, 52d31'23.70\"N)\n",
      "Band 1 Block=768x2 Type=Float32, ColorInterp=Gray\n"
     ]
    }
   ],
   "source": [
    "!gdalinfo RAD_NL25_PCP_NA_201910281900_rd_rio.tif"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "proj",
   "language": "python",
   "name": "proj"
  },
  "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.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}