Coverage for py_tools_ds/geo/coord

1# -*- coding: utf-8 -*-

3# py_tools_ds - A collection of geospatial data analysis tools that simplify standard

4# operations when handling geospatial raster and vector data as well as projections.

7# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)

8# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,

9# Germany (https://www.gfz-potsdam.de/)

10#

11# This software was developed within the context of the GeoMultiSens project funded

12# by the German Federal Ministry of Education and Research

13# (project grant code: 01 IS 14 010 A-C).

14#

15# Licensed under the Apache License, Version 2.0 (the "License");

16# you may not use this file except in compliance with the License.

17# You may obtain a copy of the License at

18#

19# http://www.apache.org/licenses/LICENSE-2.0

20#

21# Unless required by applicable law or agreed to in writing, software

22# distributed under the License is distributed on an "AS IS" BASIS,

23# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.

24# See the License for the specific language governing permissions and

25# limitations under the License.

27import numpy as np

28from shapely.geometry import box

29from shapely.geometry import Polygon # noqa F401 # flake8 issue

30from typing import Sequence, Tuple, Union # noqa F401 # flake8 issue

32from ..numeric.vector import find_nearest

33from .coord_calc import get_corner_coordinates

35__author__ = "Daniel Scheffler"

38def get_coord_grid(ULxy, LRxy, out_resXY):

39 # type: (Sequence[float, float], Sequence[float, float], Sequence[float, float]) -> np.ndarray

40 X_vec = np.arange(ULxy[0], LRxy[0], out_resXY[0])

41 Y_vec = np.arange(ULxy[1], LRxy[1], out_resXY[1])

43 # noinspection PyTypeChecker

44 return np.meshgrid(X_vec, Y_vec)

47def snap_bounds_to_pixGrid(bounds, gt, roundAlg='auto'):

48 # type: (Sequence[float], Sequence[float], str) -> Tuple[float, float, float, float]

49 """Snap the given bounds to the given grid (defined by gt) under the use of the given round algorithm.

51 NOTE: asserts equal projections of source and target grid

53 :param bounds: <tuple, list> (xmin, ymin, xmax, ymax)

54 :param gt: <tuple, list> GDAL geotransform

55 :param roundAlg: <str> 'auto', 'on', 'off'

56 :return:

57 """

58 in_xmin, in_ymin, in_xmax, in_ymax = bounds

59 xgrid = np.arange(gt[0], gt[0] + 2 * gt[1], gt[1])

60 ygrid = np.arange(gt[3], gt[3] + 2 * abs(gt[5]), abs(gt[5]))

61 xmin = find_nearest(xgrid, in_xmin, roundAlg, extrapolate=True)

62 ymax = find_nearest(ygrid, in_ymax, roundAlg, extrapolate=True)

63 xmax = find_nearest(xgrid, in_xmax, roundAlg, extrapolate=True)

64 ymin = find_nearest(ygrid, in_ymin, roundAlg, extrapolate=True)

66 return xmin, ymin, xmax, ymax

69def is_coord_grid_equal(gt, xgrid, ygrid, tolerance=0.):

70 # type: (Sequence[float], Sequence[float], Sequence[float], float) -> bool

71 """Check if a given GeoTransform exactly matches the given X/Y grid.

73 :param gt: GDAL GeoTransform

74 :param xgrid: numpy array defining the coordinate grid in x-direction

75 :param ygrid: numpy array defining the coordinate grid in y-direction

76 :param tolerance: float value defining a tolerance, e.g. 1e-8

77 :return:

78 """

79 ULx, xgsd, holder, ULy, holder, ygsd = gt

80 if all([is_point_on_grid((ULx, ULy), xgrid, ygrid, tolerance),

81 abs(xgsd - abs(xgrid[1] - xgrid[0])) <= tolerance,

82 abs(abs(ygsd) - abs(ygrid[1] - ygrid[0])) <= tolerance]):

83 return True

84 else:

85 return False

88def is_point_on_grid(pointXY, xgrid, ygrid, tolerance=0.):

89 # type: (Sequence[float, float], Sequence[float], Sequence[float], float) -> bool

90 """Check if a given point is exactly on the given coordinate grid.

92 :param pointXY: (X,Y) coordinates of the point to check

93 :param xgrid: numpy array defining the coordinate grid in x-direction

94 :param ygrid: numpy array defining the coordinate grid in y-direction

95 :param tolerance: float value defining a tolerance, e.g. 1e-8

96 """

97 if abs(find_nearest(xgrid, pointXY[0], extrapolate=True) - pointXY[0]) <= tolerance and \

98 abs(find_nearest(ygrid, pointXY[1], extrapolate=True) - pointXY[1]) <= tolerance:

99 return True

100 else:

101 return False

102

103

104def find_nearest_grid_coord(valXY, gt, rows, cols, direction='NW', extrapolate=True):

105 # type: (Sequence[float, float], Sequence[float], int, int, str, bool) -> Tuple[float, float]

106 UL, LL, LR, UR = get_corner_coordinates(gt=gt, rows=rows, cols=cols) # (x,y) tuples

107 round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SE': 'on'}[direction]

108 round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[direction]

109 tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])

110 tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))

111 tgt_x = find_nearest(tgt_xgrid, valXY[0], roundAlg=round_x, extrapolate=extrapolate)

112 tgt_y = find_nearest(tgt_ygrid, valXY[1], roundAlg=round_y, extrapolate=extrapolate)

113

114 return tgt_x, tgt_y

115

116

117def move_shapelyPoly_to_image_grid(shapelyPoly, gt, rows, cols, moving_dir='NW'):

118 # type: (Union[Polygon, box], Sequence[float], int, int, str) -> box

119 polyULxy = (min(shapelyPoly.exterior.coords.xy[0]), max(shapelyPoly.exterior.coords.xy[1]))

120 polyLRxy = (max(shapelyPoly.exterior.coords.xy[0]), min(shapelyPoly.exterior.coords.xy[1]))

121 UL, LL, LR, UR = get_corner_coordinates(gt=gt, rows=rows, cols=cols) # (x,y)

122 round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SE': 'on'}[moving_dir]

123 round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[moving_dir]

124 tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])

125 tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))

126 tgt_xmin = find_nearest(tgt_xgrid, polyULxy[0], roundAlg=round_x, extrapolate=True)

127 tgt_xmax = find_nearest(tgt_xgrid, polyLRxy[0], roundAlg=round_x, extrapolate=True)

128 tgt_ymin = find_nearest(tgt_ygrid, polyLRxy[1], roundAlg=round_y, extrapolate=True)

129 tgt_ymax = find_nearest(tgt_ygrid, polyULxy[1], roundAlg=round_y, extrapolate=True)

130

131 return box(tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax)

Coverage for py_tools_ds/geo/coord_grid.py: 49%

51 statements