coord_grid.py 6.04 KB
Newer Older
1
2
# -*- coding: utf-8 -*-

3
4
# py_tools_ds - A collection of geospatial data analysis tools that simplify standard
# operations when handling geospatial raster and vector data as well as projections.
5
#
6
7
8
9
# Copyright (C) 2016-2021
# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
#   Germany (https://www.gfz-potsdam.de/)
10
11
12
13
14
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
15
16
17
# 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
18
#
19
#   http://www.apache.org/licenses/LICENSE-2.0
20
#
21
22
23
24
25
# 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.
26

27
28
import numpy as np
from shapely.geometry import box
29
30
from shapely.geometry import Polygon  # noqa F401  # flake8 issue
from typing import Sequence, Tuple, Union  # noqa F401  # flake8 issue
31
32

from ..numeric.vector import find_nearest
33
34
35
from .coord_calc import get_corner_coordinates

__author__ = "Daniel Scheffler"
36
37


38
def get_coord_grid(ULxy, LRxy, out_resXY):
39
    # type: (Sequence[float, float], Sequence[float, float], Sequence[float, float]) -> np.ndarray
40
41
    X_vec = np.arange(ULxy[0], LRxy[0], out_resXY[0])
    Y_vec = np.arange(ULxy[1], LRxy[1], out_resXY[1])
42
43

    # noinspection PyTypeChecker
44
45
46
    return np.meshgrid(X_vec, Y_vec)


47
def snap_bounds_to_pixGrid(bounds, gt, roundAlg='auto'):
48
    # type: (Sequence[float], Sequence[float], str) -> Tuple[float, float, float, float]
Daniel Scheffler's avatar
Daniel Scheffler committed
49
50
    """Snap the given bounds to the given grid (defined by gt) under the use of the given round algorithm.

51
52
53
54
55
56
57
    NOTE: asserts equal projections of source and target grid

    :param bounds:      <tuple, list> (xmin, ymin, xmax, ymax)
    :param gt:          <tuple, list> GDAL geotransform
    :param roundAlg:    <str> 'auto', 'on', 'off'
    :return:
    """
58
59
60
    in_xmin, in_ymin, in_xmax, in_ymax = bounds
    xgrid = np.arange(gt[0], gt[0] + 2 * gt[1], gt[1])
    ygrid = np.arange(gt[3], gt[3] + 2 * abs(gt[5]), abs(gt[5]))
61
62
63
64
    xmin = find_nearest(xgrid, in_xmin, roundAlg, extrapolate=True)
    ymax = find_nearest(ygrid, in_ymax, roundAlg, extrapolate=True)
    xmax = find_nearest(xgrid, in_xmax, roundAlg, extrapolate=True)
    ymin = find_nearest(ygrid, in_ymin, roundAlg, extrapolate=True)
65

66
67
68
    return xmin, ymin, xmax, ymax


69
def is_coord_grid_equal(gt, xgrid, ygrid, tolerance=0.):
70
    # type: (Sequence[float], Sequence[float], Sequence[float], float) -> bool
Daniel Scheffler's avatar
Daniel Scheffler committed
71
    """Check if a given GeoTransform exactly matches the given X/Y grid.
Daniel Scheffler's avatar
Daniel Scheffler committed
72
73
74
75
76
77
78

    :param gt:          GDAL GeoTransform
    :param xgrid:       numpy array defining the coordinate grid in x-direction
    :param ygrid:       numpy array defining the coordinate grid in y-direction
    :param tolerance:   float value defining a tolerance, e.g. 1e-8
    :return:
    """
79
80
81
82
    ULx, xgsd, holder, ULy, holder, ygsd = gt
    if all([is_point_on_grid((ULx, ULy), xgrid, ygrid, tolerance),
            abs(xgsd - abs(xgrid[1] - xgrid[0])) <= tolerance,
            abs(abs(ygsd) - abs(ygrid[1] - ygrid[0])) <= tolerance]):
83
84
85
86
87
        return True
    else:
        return False


Daniel Scheffler's avatar
Daniel Scheffler committed
88
def is_point_on_grid(pointXY, xgrid, ygrid, tolerance=0.):
89
    # type: (Sequence[float, float], Sequence[float], Sequence[float], float) -> bool
Daniel Scheffler's avatar
Daniel Scheffler committed
90
    """Check if a given point is exactly on the given coordinate grid.
Daniel Scheffler's avatar
Daniel Scheffler committed
91

92
93
94
    :param pointXY:     (X,Y) coordinates of the point to check
    :param xgrid:       numpy array defining the coordinate grid in x-direction
    :param ygrid:       numpy array defining the coordinate grid in y-direction
Daniel Scheffler's avatar
Daniel Scheffler committed
95
    :param tolerance:   float value defining a tolerance, e.g. 1e-8
96
    """
97
98
    if abs(find_nearest(xgrid, pointXY[0], extrapolate=True) - pointXY[0]) <= tolerance and \
            abs(find_nearest(ygrid, pointXY[1], extrapolate=True) - pointXY[1]) <= tolerance:
99
100
101
102
103
        return True
    else:
        return False


104
def 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
108
    round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SE': 'on'}[direction]
    round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[direction]
109
110
111
112
    tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])
    tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))
    tgt_x = find_nearest(tgt_xgrid, valXY[0], roundAlg=round_x, extrapolate=extrapolate)
    tgt_y = find_nearest(tgt_ygrid, valXY[1], roundAlg=round_y, extrapolate=extrapolate)
113

114
    return tgt_x, tgt_y
115
116


117
118
def move_shapelyPoly_to_image_grid(shapelyPoly, gt, rows, cols, moving_dir='NW'):
    # type: (Union[Polygon, box], Sequence[float], int, int, str) -> box
119
120
    polyULxy = (min(shapelyPoly.exterior.coords.xy[0]), max(shapelyPoly.exterior.coords.xy[1]))
    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
123
124
125
    round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SE': 'on'}[moving_dir]
    round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[moving_dir]
    tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])
    tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))
126
127
128
129
    tgt_xmin = find_nearest(tgt_xgrid, polyULxy[0], roundAlg=round_x, extrapolate=True)
    tgt_xmax = find_nearest(tgt_xgrid, polyLRxy[0], roundAlg=round_x, extrapolate=True)
    tgt_ymin = find_nearest(tgt_ygrid, polyLRxy[1], roundAlg=round_y, extrapolate=True)
    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)