Source code for nsaph_gis.geometry

"""
Utilities to work with points in a raster
"""

#  Copyright (c) 2021. Harvard University
#
#  Developed by Research Software Engineering,
#  Faculty of Arts and Sciences, Research Computing (FAS RC)
#  Author: Michael A Bouzinier
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#         http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#

import itertools
from typing import Optional

from numpy.ma import masked
from rasterstats import point


[docs]class PointInRaster: """ Class denoting a point in a raster. This a wrapper around class rasterstats.point optimizing some operations, primarily the bilinear interpolation See also https://pythonhosted.org/rasterstats/_modules/rasterstats/point.html """ COMPLETELY_MASKED = 1 PARTIALLY_MASKED = 2 def __init__(self, raster, affine, x, y): self.x = None '''X coordinate of the point''' self.y = None '''Y coordinate of the point''' self.window = None '''Window representing 2x2 window whose center points encompass point''' self.window, unitxy = point.point_window_unitxy(x, y, affine) self.x, self.y = unitxy self.masked = 0 m = 0 array = raster.read(window=self.window, masked=True).array for i, j in itertools.product([0,1], [0,1]): r = self.window[0][0] + i c = self.window[1][0] + j if array[i, j] is masked: m += 1 elif raster.array[r, c] is masked: m += 1 else: self.r, self.c = r, c if m == 4: self.masked = self.COMPLETELY_MASKED elif m > 0: self.masked = self.PARTIALLY_MASKED
[docs] def is_masked(self) -> bool: return self.masked == self.COMPLETELY_MASKED
[docs] def array(self, raster): return raster.array[ self.window[0][0]:self.window[0][1], self.window[1][0]:self.window[1][1], ]
[docs] def bilinear(self, raster) -> Optional[float]: """ An optimized version of rasterstats.point function: given a point's window as 2x2 array, and x, y as its coordinates, treat center points as a unit square. +---+---+ | A | B | +----+ +---+---+ => | | | C | D | +----+ +---+---+ e.g.: Center of A is at (0, 1) on unit square, D is at (1, 0), etc :param raster: Raster, to which the point belongs :return: the value for the fractional row/col using bilinear interpolation between the cells """ if self.masked == self.COMPLETELY_MASKED: return None if self.masked == self.PARTIALLY_MASKED: return raster.array[self.r, self.c] array = self.array(raster) x, y = self.x, self.y ulv, urv, llv, lrv = array[0,0], array[0,1], array[1,0], array[1,1] return ( (llv * (1 - x) * (1 - y)) + (lrv * x * (1 - y)) + (ulv * (1 - x) * y) + (urv * x * y) )