forked from nasaharvest/kenya-crop-mask
-
Notifications
You must be signed in to change notification settings - Fork 0
/
split_tiff.py
75 lines (57 loc) · 2.47 KB
/
split_tiff.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
"""
https://gis.stackexchange.com/questions/306861/split-geotiff-into-multiple-cells-with-
"""
from shapely import geometry
from pathlib import Path
import math
import rasterio
from rasterio.mask import mask
# Takes a dataset and splits it into squares of dimensions squareDim * squareDim
def splitImageIntoCells(img_path: Path, filename: str, squareDim: int, output_folder: Path):
img = rasterio.open(img_path)
numberOfCellsWide = math.ceil(img.shape[1] / squareDim)
numberOfCellsHigh = math.ceil(img.shape[0] / squareDim)
count = 0
for hc in range(numberOfCellsHigh):
y = min(hc * squareDim, img.shape[0])
for wc in range(numberOfCellsWide):
x = min(wc * squareDim, img.shape[1])
geom = getTileGeom(img.transform, x, y, squareDim)
getCellFromGeom(img, geom, filename, count, output_folder)
count = count + 1
# Generate a bounding box from the pixel-wise coordinates using the original datasets transform property
def getTileGeom(transform, x, y, squareDim):
corner1 = (x, y) * transform
corner2 = (x + squareDim, y + squareDim) * transform
return geometry.box(corner1[0], corner1[1], corner2[0], corner2[1])
# Crop the dataset using the generated box and write it out as a GeoTIFF
def getCellFromGeom(img, geom, filename, count, output_folder):
crop, cropTransform = mask(img, [geom], crop=True)
writeImageAsGeoTIFF(
crop, cropTransform, img.meta, img.crs, f"{count}-{filename}", output_folder
)
# Write the passed in dataset as a GeoTIFF
def writeImageAsGeoTIFF(img, transform, metadata, crs, filename, output_folder):
metadata.update(
{
"driver": "GTiff",
"height": img.shape[1],
"width": img.shape[2],
"transform": transform,
"crs": crs,
}
)
with rasterio.open(output_folder / f"{filename}.tif", "w", **metadata) as dest:
dest.write(img)
if __name__ == "__main__":
images = Path("data//test").glob("*.tif")
output_folder = Path("data//test")
for idx, image in enumerate(images):
print(f"Splitting {image}")
name, start_date, end_bit = image.name.split("_")
end_date = end_bit[:10]
tile_identifier = end_bit[10:-4]
new_filename = f"{idx}-{name}{tile_identifier}_{start_date}_{end_date}"
splitImageIntoCells(image, new_filename, 1000, output_folder)
print(f"Finished {image}. Removing the original file")
image.unlink()