-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathpoints_to_raster.py
More file actions
44 lines (35 loc) · 2.1 KB
/
points_to_raster.py
File metadata and controls
44 lines (35 loc) · 2.1 KB
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
"""
Creates a raster from a point ShapeFile.
This is used in the campsite analysis to create a raster from a ShapeFile that
is generated by merging the vertices of the campsite polyline vertices and points
from a single corgrid text file.
"""
from typing import Tuple
import os
from subprocess import call, PIPE
from raster import delete_raster
from logger import Logger
def points_to_raster(gdal_grid_path: str, points_shapefile: str, z_field: str, out_raster: str, cell_size: float, buffered_extent: Tuple[float, float, float, float]) -> None:
"""
:param gdal_grid_path: The path to the GDAL Grid executable
:param points_raster: The path to the input points ShapeFile
:param out_raster: The path to the output raster
:param clip_shape_file: The path to the polygon ShapeFile to use for clipping
:param where_clause: Feature filter for selecting which features to use for clipping
https://gdal.org/programs/gdal_grid.html
"""
log = Logger('Points to Raster')
assert os.path.isfile(gdal_grid_path), f'Missing GDAL Grid executable at {gdal_grid_path}'
assert os.path.isfile(points_shapefile), f'Missing point ShapeFile operation input at {points_shapefile}'
# Make sure the rasters get removed before they get re-made
delete_raster(out_raster)
# Reset the where parameter to an empty string if no where clause is provided
# TODO: This is giving us 64-bit rasters for some reason and a weird nodata value with nan as well. We're probably losing precision somewhere
x_extent = f'{buffered_extent[0]} {buffered_extent[1]}'
y_extent = f'{buffered_extent[2]} {buffered_extent[3]}'
gdal_args = f' -a linear -zfield {z_field} -tr {cell_size} {cell_size} -txe {x_extent} -tye {y_extent} {points_shapefile} {out_raster}'
log.debug('RUNNING GdalGrid: ' + gdal_grid_path + gdal_args)
if ' ' in gdal_grid_path:
gdal_grid_path = f'"{gdal_grid_path}"'
return_val = call(gdal_grid_path + gdal_args, stdout=PIPE, shell=True)
assert return_val == 0, f'Error creating raster from points. Input points ShapeFile {points_shapefile} and output raster {out_raster}'