make_small_geotiff_landsat.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CustomArray(xr.DataArray):\n",
    "    def to_rasterio(self, filename, parse_transform=True):\n",
    "        \"\"\"\n",
    "        Write a DataArray to a file using rasterio. DataArray must have \n",
    "        dimensions x, y and band.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        filename : str\n",
    "            Path to the file to open.\n",
    "        parse_transform : bool, optional\n",
    "            Whether to parse the ``transform`` attribute out of the x and y \n",
    "            coordinates or not. The default is to automatically\n",
    "            parse the coordinates only if they are rectilinear (1D).\n",
    "            It can be useful to set ``parse_transform=False``\n",
    "            if your files are very large or if you don't need the transform.\n",
    "        \"\"\"\n",
    "        import rasterio\n",
    "        \n",
    "        if self.ndim == 2:\n",
    "            nband = 1\n",
    "            band_indicies = 1\n",
    "        else:\n",
    "            nband = self.band.shape[0]\n",
    "            band_indicies = list(range(1, nband + 1))\n",
    "        \n",
    "        ny = self.y.shape[0]\n",
    "        nx = self.x.shape[0]\n",
    "    \n",
    "        if parse_transform:\n",
    "            x_res = (self.x[1] - self.x[0]).item()\n",
    "            y_res = (self.y[1] - self.y[0]).item()\n",
    "\n",
    "            x_offset = self.x[0].item() - x_res/2.\n",
    "            y_offset = self.y[0].item() - y_res/2.\n",
    "\n",
    "            transform = rasterio.Affine(x_res, 0, x_offset, 0, y_res, y_offset)\n",
    "        else:\n",
    "            transform = None\n",
    "\n",
    "        with rasterio.open(filename, 'w+',\n",
    "                           driver='GTiff',\n",
    "                           height=ny, width=nx,\n",
    "                           dtype=str(self.dtype), count=nband,\n",
    "                           transform=transform, crs=self.crs) as dst:\n",
    "            dst.write(self.values, band_indicies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = 'landsat/data/MERCATOR_LC80210392016114LGN00_B%s.TIF'\n",
    "new_path = 'test_data/landsat/MERCATOR_LC80210392016114LGN00_B%s.TIF'\n",
    "bands = list(range(1, 12)) + ['QA']\n",
    "bands = [xr.open_rasterio(file_path%band).load()[0] for band in bands]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for band in range(1, 12):\n",
    "    original = bands[band-1]\n",
    "    new = original[3500:5000:10, 3000:4500:10]\n",
    "    new = CustomArray(new)\n",
    "    new.to_rasterio(new_path%band)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band = 'QA'\n",
    "original = bands[11]\n",
    "new = original[3500:5000:10, 3000:4500:10]\n",
    "new = CustomArray(new)\n",
    "new.to_rasterio(new_path%band)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}