Coverage for arosics/geometry.py: 89%

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

2 

3# AROSICS - Automated and Robust Open-Source Image Co-Registration Software 

4# 

5# Copyright (C) 2017-2024 

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

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

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

9# 

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

11# by the German Federal Ministry of Education and Research 

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

13# 

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

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

16# You may obtain a copy of the License at 

17# 

18# https://www.apache.org/licenses/LICENSE-2.0 

19# 

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

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

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

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

24# limitations under the License. 

25 

26import warnings 

27import sys 

28import os 

29from typing import Union 

30 

31# custom 

32import numpy as np 

33from geopandas import GeoDataFrame 

34 

35# internal modules 

36from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions 

37from py_tools_ds.geo.coord_trafo import pixelToMapYX, imYX2mapYX 

38from py_tools_ds.geo.raster.reproject import warp_ndarray 

39from geoarray import GeoArray 

40 

41__author__ = 'Daniel Scheffler' 

42 

43 

44def angle_to_north(XY): 

45 """Calculate the angle in degrees of a given line to north in clockwise direction. 

46 

47 Angle definition: between [origin:[0,0],north:[0,1]] and [origin:[0,0],pointXY:[X,Y]]. 

48 """ 

49 XY = np.array(XY) 

50 XYarr = XY if len(XY.shape) == 2 else XY.reshape((1, 2)) 

51 return np.abs(np.degrees(np.array(np.arctan2(XYarr[:, 1], XYarr[:, 0]) - np.pi / 2)) % 360) 

52 

53 

54def get_true_corner_mapXY(fPath_or_geoarray, band=0, noDataVal=None, mp=1, v=0, q=0): # pragma: no cover 

55 """Return the actual map corner coordinates of a given image file or GeoArray instance. 

56 

57 :param fPath_or_geoarray: 

58 :param band: <int> index of the band to be used (starting with 0) 

59 :param noDataVal: 

60 :param mp: 

61 :param v: 

62 :param q: 

63 :return: 

64 """ 

65 # FIXME this function is not used anymore 

66 warnings.warn('This function is not in use anymore. Use it on your own risk!', DeprecationWarning) 

67 geoArr = GeoArray(fPath_or_geoarray) if not isinstance(fPath_or_geoarray, GeoArray) else fPath_or_geoarray 

68 

69 rows, cols = geoArr.shape[:2] 

70 gt, prj = geoArr.geotransform, geoArr.projection 

71 

72 assert gt and prj, 'GeoTransform an projection must be given for calculation of LonLat corner coordinates.' 

73 

74 mask_1bit = np.zeros((rows, cols), dtype='uint8') # zeros -> image area later overwritten by ones 

75 

76 if noDataVal is None: 

77 mask_1bit[:, :] = 1 

78 elif noDataVal == 'ambiguous': 

79 warnings.warn("No data value could not be automatically detected. Thus the matching window used for shift " 

80 "calculation had to be centered in the middle of the overlap area without respecting no data " 

81 "values. To avoid this provide the correct no data values for reference and shift image via " 

82 "'-nodata'") 

83 mask_1bit[:, :] = 1 

84 else: 

85 band_data = geoArr[band] # TODO implement gdal_ReadAsArray_mp (reading in multiprocessing) 

86 mask_1bit[band_data != noDataVal] = 1 

87 

88 if v: 

89 print('detected no data value', noDataVal) 

90 

91 try: 

92 corner_coords_YX = calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='shapely') 

93 except Exception: 

94 if v: 

95 warnings.warn("\nCalculation of corner coordinates failed within algorithm 'shapely' (Exception: %s)." 

96 " Using algorithm 'numpy' instead." % sys.exc_info()[1]) 

97 # FIXME numpy algorithm returns wrong values for S2A_OPER_MSI_L1C_TL_SGS__20160608T153121_A005024_T33UUU_B03.jp2 

98 # FIXME (Hannes) 

99 corner_coords_YX = \ 

100 calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='numpy') 

101 

102 if len(corner_coords_YX) == 4: # this avoids shapely self intersection 

103 corner_coords_YX = list(np.array(corner_coords_YX)[[0, 1, 3, 2]]) # UL, UR, LL, LR => UL, UR, LR, LL 

104 

105 # check if enough unique coordinates have been found 

106 if not len(GeoDataFrame(corner_coords_YX).drop_duplicates().values) >= 3: 

107 if not q: 

108 warnings.warn('\nThe algorithm for automatically detecting the actual image coordinates did not find ' 

109 'enough unique corners. Using outer image corner coordinates instead.') 

110 corner_coords_YX = ((0, 0), (0, cols - 1), (rows - 1, 0), (rows - 1, cols - 1)) 

111 

112 # check if all points are unique 

113 # all_coords_are_unique = len([UL, UR, LL, LR]) == len(GeoDataFrame([UL, UR, LL, LR]).drop_duplicates().values) 

114 # UL, UR, LL, LR = \ 

115 # (UL, UR, LL, LR) if all_coords_are_unique else ((0, 0), (0, cols-1), (rows-1, 0), (rows-1, cols-1)) 

116 

117 def get_mapYX(YX): return pixelToMapYX(list(reversed(YX)), geotransform=gt, projection=prj)[0] 

118 corner_pos_XY = [list(reversed(i)) for i in [get_mapYX(YX) for YX in corner_coords_YX]] 

119 return corner_pos_XY 

120 

121 

122def get_subset_GeoTransform(gt_fullArr, subset_box_imYX): 

123 gt_subset = list(gt_fullArr[:]) # copy 

124 gt_subset[3], gt_subset[0] = imYX2mapYX(subset_box_imYX[0], gt_fullArr) 

125 return gt_subset 

126 

127 

128def get_gdalReadInputs_from_boxImYX(boxImYX): 

129 """Return row_start,col_start,rows_count,cols_count and assumes boxImYX as [UL_YX,UR_YX,LR_YX,LL_YX).""" 

130 rS, cS = boxImYX[0] 

131 clip_sz_x = abs(boxImYX[1][1] - boxImYX[0][1]) # URx-ULx 

132 clip_sz_y = abs(boxImYX[0][0] - boxImYX[3][0]) # ULy-LLy 

133 return cS, rS, clip_sz_x, clip_sz_y 

134 

135 

136def get_GeoArrayPosition_from_boxImYX(boxImYX): 

137 """Return row_start,row_end,col_start,col_end and assumes boxImYX as [UL_YX,UR_YX,LR_YX,LL_YX).""" 

138 rS, cS = boxImYX[0] # UL 

139 rE, cE = boxImYX[2] # LR 

140 return rS, rE - 1, cS, cE - 1 # -1 because boxImYX represents outer box and includes the LR corner of LR pixel 

141 

142 

143def has_metaRotation(path_or_geoarray: Union[GeoArray, str]): 

144 """Return True if there is a row or column rotation due to the given GDAL GeoTransform tuple.""" 

145 gt = GeoArray(path_or_geoarray).gt 

146 

147 return gt[2] or gt[4] 

148 

149 

150def remove_metaRotation(gA_rot: GeoArray, rspAlg='cubic') -> GeoArray: 

151 """Remove any metadata rotation (a rotation that only exists in the map info).""" 

152 gA = GeoArray(*warp_ndarray(gA_rot[:], gA_rot.gt, gA_rot.prj, 

153 rspAlg=rspAlg, 

154 # out_gsd=(gA_rot.xgsd, gA_rot.ygsd) 

155 ), 

156 nodata=gA_rot.nodata) 

157 gA.basename = os.path.basename(gA.basename) 

158 gA.meta = gA.meta 

159 

160 return gA