makegeotiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "16919fb0",
   "metadata": {},
   "source": [
    "# reference\n",
    "\n",
    " - https://club.informatix.co.jp/?p=1293"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdce0fee",
   "metadata": {},
   "source": [
    "# environment\n",
    "\n",
    "```sh\n",
    "conda create -y -n tmp python=3.7\n",
    "conda activate tmp \n",
    "conda install -y -c conda-forge rioxarray\n",
    "conda install -y -c conda-forge gdal\n",
    "conda install -y -c conda-forge jupyter\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ae24a97",
   "metadata": {},
   "source": [
    "# make geotiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d4c7a44f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "import glob \n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import rioxarray\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5a3e3502",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getcoords(word, xorg, yorg, dx, dy):\n",
    "    anum = lambda x : ord(x) - ord('a')\n",
    "    \n",
    "    v = str(word)\n",
    "    \n",
    "    if re.match('\\d+', v):\n",
    "        iy, ix = int(v[0]), int(v[1])\n",
    "    else:\n",
    "        iy, ix = anum(v[0]), anum(v[1])\n",
    "        \n",
    "    yc = yorg - iy*dy\n",
    "    xc = xorg + ix*dx\n",
    "    \n",
    "    return xc, yc\n",
    "\n",
    "# 4分割\n",
    "def getcoords2(word, xorg, yorg, dx, dy):\n",
    "    v = int(word)\n",
    "    \n",
    "    if v == 1:\n",
    "        iy, ix = 0, 0\n",
    "    elif v == 2:\n",
    "        iy, ix = 0, 1\n",
    "    elif v == 3:\n",
    "        iy, ix = 1, 0\n",
    "    elif v == 4:\n",
    "        iy, ix = 1, 1\n",
    "    else:\n",
    "        print('error')\n",
    "        \n",
    "    yc = yorg - iy*dy\n",
    "    xc = xorg + ix*dx\n",
    "        \n",
    "    return xc, yc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "76ad70ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs  =glob.glob('*.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "847ba53f",
   "metadata": {},
   "outputs": [],
   "source": [
    "for fp in fs:\n",
    "\n",
    "    f = os.path.basename(fp)\n",
    "    \n",
    "    df = pd.read_csv(fp, header=None)\n",
    "    X = df[1].values\n",
    "    Y = df[2].values\n",
    "    Z = df[3].values\n",
    "    \n",
    "   # float to int \n",
    "    Z = (100*Z).astype(np.int32)\n",
    "    \n",
    "    w = f[:8]\n",
    "    nepsg = int(w[:2])\n",
    "    \n",
    "    # 図郭の北西の座標を求める\n",
    "    # 地図情報レベル 50000\n",
    "    xc, yc = getcoords(w[2:4], xorg=-160000, yorg=300000, dx=40000, dy=30000)\n",
    "    # 地図情報レベル 5000\n",
    "    xc, yc = getcoords(w[4:6], xorg=xc, yorg=yc, dx=4000, dy=3000)\n",
    "    # 地図情報レベル 5000を4分割\n",
    "    xc, yc = getcoords2(w[6], xorg=xc, yorg=yc, dx=2000, dy=1500)\n",
    "#     # # さらに4分割\n",
    "#     xc, yc = getcoords2(w[7], xorg=xc, yorg=yc, dx=1000, dy=750)\n",
    "    \n",
    "    delta = float(0.5)\n",
    "    dx=2000\n",
    "    dy=1500\n",
    "    \n",
    "    ix = ((X - xc - 0.5*delta)/delta).astype(int)\n",
    "    iy = ((- Y + yc - 0.5*delta)/delta).astype(int)\n",
    "    \n",
    "    zout= np.full((int(dx/delta),int(dy/delta)), int(-9999) )\n",
    "#     zout= np.full((int(dx/delta),int(dy/delta)), float(-9999) )\n",
    "    \n",
    "    for ixp, iyp, zp in zip(ix, iy, Z) : zout[ixp,iyp] = zp\n",
    "    \n",
    "    xcoord = np.arange(xc, xc+dx, delta) + 0.5*delta\n",
    "    ycoord = np.arange(yc, yc-dy, -delta) - 0.5*delta\n",
    "    \n",
    "    epsg = str(6668 + nepsg)\n",
    "    \n",
    "    ds = xr.Dataset({'z': (['y','x'], zout.T) }, coords={'x': xcoord, 'y': ycoord}) \n",
    "    ds = ds.rio.write_crs('EPSG:' + epsg, inplace = True)\n",
    "    # ds.rio.reproject(\"epsg:****\")\n",
    "    \n",
    "    # export geotiff file\n",
    "    out = ds['z'].rio.to_raster(  f[:-4] + '.tif')\n",
    "    del out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "def33072",
   "metadata": {},
   "source": [
    "# make vrt file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "580e3d24",
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal \n",
    "\n",
    "opt=gdal.BuildVRTOptions(VRTNodata=float(-9999), srcNodata=float(-9999))\n",
    "my_vrt = gdal.BuildVRT('output.vrt', glob.glob( '*.tif'), options=opt)\n",
    "del my_vrt"
   ]
  }
 ],
 "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.7.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}