Skip to content

Instantly share code, notes, and snippets.

@neon-ninja

neon-ninja/get_consistent_extents.ipynb

Created January 5, 2022 03:11
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save neon-ninja/73ab2c84541a9bc8e7a67a96805db449 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save neon-ninja/73ab2c84541a9bc8e7a67a96805db449 to your computer and use it in GitHub Desktop.
Download ZIP
Find the common extent among GeoTIFFs, clip them all to the common extent, reproject to common transform
Display the source blob
Display the rendered blob
Raw
get_consistent_extents.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "6ee3da8d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['Epoch1.tif', 'Epoch10.tif', 'Epoch11.tif', 'Epoch2.tif', 'Epoch3.tif', 'Epoch4v2.tif', 'Epoch5.tif', 'Epoch6.tif', 'Epoch7.tif', 'Epoch9.tif', 'New_Epoch8.tif']\n"
]
}
],
"source": [
"import rasterio as rio\n",
"from glob import glob\n",
"import pandas as pd\n",
"from rasterio.mask import mask\n",
"from rasterio.warp import reproject, calculate_default_transform\n",
"from shapely.geometry import box\n",
"import shapely\n",
"import os\n",
"from tqdm.auto import tqdm\n",
"files = sorted(f for f in glob(\"*.tif\") if \"bgis\" not in f)\n",
"print(files)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "961f49fa",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[BoundingBox(left=1654871.0, bottom=5315335.0, right=1658478.0, top=5323188.0),\n",
" BoundingBox(left=1654870.5304000105, bottom=5315334.320583384, right=1658477.5304000105, top=5323187.320583384),\n",
" BoundingBox(left=1654871.0, bottom=5315335.0, right=1658478.0, top=5323188.0),\n",
" BoundingBox(left=1654871.08, bottom=5315334.734, right=1658478.08, top=5323187.734),\n",
" BoundingBox(left=1654871.0, bottom=5315335.0, right=1658478.0, top=5323188.0),\n",
" BoundingBox(left=1654871.0, bottom=5315335.0, right=1658478.0, top=5323188.0),\n",
" BoundingBox(left=1654872.0, bottom=5315336.0, right=1658477.0, top=5323187.0),\n",
" BoundingBox(left=1654870.98, bottom=5315335.0, right=1658477.98, top=5323188.0),\n",
" BoundingBox(left=1654871.316156081, bottom=5315335.087024724, right=1658477.316156081, top=5323188.087024724),\n",
" BoundingBox(left=1654870.4495669291, bottom=5315334.49115541, right=1658477.4495669291, top=5323187.49115541),\n",
" BoundingBox(left=1654870.5350000001, bottom=5315335.012, right=1658477.5350000001, top=5323188.012)]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tiles = [rio.open(f) for f in files]\n",
"bounds = [f.bounds for f in tiles]\n",
"bounds"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "9d7cd56a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1654872.0, 5315336.0, 1658477.0, 5323187.0]\n"
]
},
{
"data": {
"text/plain": [
"{'type': 'Polygon',\n",
" 'coordinates': (((1658477.0, 5315336.0),\n",
" (1658477.0, 5323187.0),\n",
" (1654872.0, 5323187.0),\n",
" (1654872.0, 5315336.0),\n",
" (1658477.0, 5315336.0)),)}"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df=pd.DataFrame(bounds)\n",
"common_extent = [df.left.max(), df.bottom.max(), df.right.min(), df.top.min()]\n",
"print(common_extent)\n",
"bbox = box(*common_extent)\n",
"bbox = shapely.geometry.mapping(bbox)\n",
"bbox"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8aed687c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028230607370965e+38, 'width': 3607, 'height': 7853, 'count': 1, 'transform': Affine(1.0, 0.0, 1654871.0,\n",
" 0.0, -1.0, 5323188.0)}\n",
"{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028230607370965e+38, 'width': 3605, 'height': 7851, 'count': 1, 'transform': Affine(1.0, 0.0, 1654872.0,\n",
" 0.0, -1.0, 5323187.0)}\n"
]
}
],
"source": [
"src = tiles[0]\n",
"ref_img, ref_transform = mask(src, shapes=[bbox], crop=True)\n",
"ref_img.shape, ref_transform\n",
"print({k: v for k,v in src.meta.items() if k!=\"crs\"})\n",
"out_meta = src.meta.copy()\n",
"out_meta.update({\n",
" \"height\": ref_img.shape[1],\n",
" \"width\": ref_img.shape[2],\n",
" \"transform\": ref_transform\n",
"})\n",
"print({k: v for k,v in out_meta.items() if k!=\"crs\"})"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "a79f532e",
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "f5a591ae58de4eb19ca431396e4c2cef",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/11 [00:00<?, ?it/s]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Epoch1.tif\n",
"Epoch10.tif\n",
"Epoch11.tif\n",
"Epoch2.tif\n",
"Epoch3.tif\n",
"Epoch4v2.tif\n",
"Epoch5.tif\n",
"Epoch6.tif\n",
"Epoch7.tif\n",
"Epoch9.tif\n",
"New_Epoch8.tif\n"
]
}
],
"source": [
"os.makedirs(\"clip\", exist_ok=True)\n",
"for f in tqdm(files):\n",
" print(f)\n",
" with rio.open(f) as src:\n",
" # Clip to bbox\n",
" out_img, out_transform = mask(src, shapes=[bbox], crop=True)\n",
" with rio.open(f\"clip/{f}\", \"w\", **out_meta) as dst:\n",
" # Reproject clip result to reference transform\n",
" reproject(\n",
" source=out_img,\n",
" destination=rio.band(dst, 1),\n",
" src_transform=out_transform,\n",
" src_crs=src.crs,\n",
" dst_transform=dst.transform,\n",
" dst_crs=dst.crs)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f7b0aefd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028230607370965e+38, 'width': 3605, 'height': 7851, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(1.0, 0.0, 1654872.0,\n",
" 0.0, -1.0, 5323187.0)}\n"
]
}
],
"source": [
"tiles = [rio.open(f\"clip/{f}\") for f in files]\n",
"print(tiles[0].meta)\n",
"for t in tiles:\n",
" assert t.meta == tiles[0].meta"
]
}
],
"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.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment