Coverage for sarvey/osm_utils.py: 0%

1#!/usr/bin/env python 

2 

3# SARvey - A multitemporal InSAR time series tool for the derivation of displacements. 

4# 

5# Copyright (C) 2021-2024 Andreas Piter (IPI Hannover, piter@ipi.uni-hannover.de) 

6# 

7# This software was developed together with FERN.Lab (fernlab@gfz-potsdam.de) in the context 

8# of the SAR4Infra project with funds of the German Federal Ministry for Digital and 

9# Transport and contributions from Landesamt fuer Vermessung und Geoinformation 

10# Schleswig-Holstein and Landesbetrieb Strassenbau und Verkehr Schleswig-Holstein. 

11# 

12# This program is free software: you can redistribute it and/or modify it under 

13# the terms of the GNU General Public License as published by the Free Software 

14# Foundation, either version 3 of the License, or (at your option) any later 

15# version. 

16# 

17# Important: This package uses PyMaxFlow. The core of PyMaxflows library is the C++ 

18# implementation by Vladimir Kolmogorov. It is also licensed under the GPL, but it REQUIRES that you 

19# cite [BOYKOV04] (see LICENSE) in any resulting publication if you use this code for research purposes. 

20# This requirement extends to SARvey. 

21# 

22# This program is distributed in the hope that it will be useful, but WITHOUT 

23# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

24# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

25# details. 

26# 

27# You should have received a copy of the GNU Lesser General Public License along 

28# with this program. If not, see <https://www.gnu.org/licenses/>. 

29 

30"""Osm utils module for SARvey.""" 

31import numpy as np 

32import overpy 

33from logging import Logger 

34from shapely import Point 

35 

36from mintpy.utils import readfile, utils as ut 

37 

38 

39def getSpatialExtend(*, geom_file: str, logger: Logger): 

40 """Get spatial extend of the radar image. 

41 

42 Parameters 

43 ---------- 

44 geom_file: str 

45 path of geometryRadar.h5 file 

46 logger: Logger 

47 Logging handler. 

48 

49 Returns 

50 ------- 

51 ll_corner_wgs: list 

52 list of coordinates of the lower-left corner of the radar image in WGS84 coordinates. 

53 ur_corner_wgs: list 

54 list of coordinates of the upper-right corner of the radar image in WGS84 coordinates. 

55 coord: np.ndarray 

56 coordinates of all pixels in the radar image in WGS84. 

57 atr: dict 

58 metadata dictionary from geometryRadar.h5. 

59 """ 

60 logger.info(msg='read spatial extend from geometryRadar.h5') 

61 _, atr = readfile.read(geom_file) 

62 coord = ut.coordinate(atr, lookup_file=geom_file) 

63 lat, atr = readfile.read(geom_file, datasetName='latitude') 

64 lon, _ = readfile.read(geom_file, datasetName='longitude') 

65 

66 # radar image is flipped up-down 

67 # unclear: check if bounding box fits. Otherwise, change to max and min values of lat and lon 

68 ll_bbox = [np.nanmin(lat), np.nanmin(lon)] 

69 ur_bbox = [np.nanmax(lat), np.nanmax(lon)] 

70 

71 img_ext = [ 

72 Point(lon[0, 0], lat[0, 0]), 

73 Point(lon[-1, 0], lat[-1, 0]), 

74 Point(lon[-1, -1], lat[-1, -1]), 

75 Point(lon[0, -1], lat[0, -1]) 

76 ] 

77 return ll_bbox, ur_bbox, img_ext, coord, atr 

78 

79 

80def runOsmQuery(*, ll_corner_wgs: np.ndarray, ur_corner_wgs: np.ndarray, type_list: list, 

81 logger: Logger) -> overpy.Result: 

82 """Query OSM database for transport infrastructure within the spatial extent of the radar image. 

83 

84 Parameters 

85 ---------- 

86 ll_corner_wgs: np.ndarray 

87 coordinates of the lower-left corner of the radar image in WGS84 coordinates. 

88 ur_corner_wgs: np.ndarray 

89 coordinates of the upper-right corner of the radar image in WGS84 coordinates. 

90 type_list: list 

91 List of street types that shall be queried at the OSM database. 

92 logger: Logger 

93 Logging handler. 

94 

95 Returns 

96 ------- 

97 result: overpy.Result 

98 results of the overpy query to OSM database. 

99 """ 

100 # Initialize overpass connection 

101 api = overpy.Overpass() 

102 

103 # Request data from API 

104 logger.info(msg='querying OSM database for infra types...') 

105 # query_cmd = "way({},{},{},{}) [""highway=motorway_link""]; (._;>;); out body;" 

106 

107 query_cmd = "[bbox: {},{},{},{}];(" 

108 for infra_type in type_list: 

109 logger.info(msg='\t - {}'.format(infra_type)) 

110 if infra_type == 'rail': 

111 query_cmd += "way[railway={}];".format(infra_type) 

112 else: 

113 query_cmd += "way[highway={}];".format(infra_type) 

114 

115 query_cmd += ");(._; >;); out body;" # skel 

116 

117 cmd = query_cmd.format(ll_corner_wgs[0], ll_corner_wgs[1], 

118 ur_corner_wgs[0], ur_corner_wgs[1]) 

119 logger.info(msg="\n" + cmd + "\n") 

120 result = api.query(cmd) 

121 

122 if len(result.ways) == 0: 

123 logger.error(msg='Empty OSM query results. No roads or railway tracks found.') 

124 raise ValueError 

125 

126 logger.info(msg='...done.') 

127 return result 

128 

129 

130def runOsmQueryBridge(*, ll_corner_wgs: np.ndarray, ur_corner_wgs: np.ndarray, bridge_highway: bool, 

131 bridge_railway: bool, logger: Logger) -> overpy.Result: 

132 """Query OSM database for bridges of transport infrastructure within the spatial extent of the radar image. 

133 

134 Parameters 

135 ---------- 

136 ll_corner_wgs: np.ndarray 

137 coordinates of the lower-left corner of the radar image in WGS84 coordinates. 

138 ur_corner_wgs: np.ndarray 

139 coordinates of the upper-right corner of the radar image in WGS84 coordinates. 

140 bridge_highway: bool 

141 Set true to query highway bridges. 

142 bridge_railway: bool 

143 Set true to query railway bridges. 

144 logger: Logger 

145 Logging handler. 

146 

147 Returns 

148 ------- 

149 result: overpy.Result 

150 results of the overpy query to OSM database. 

151 """ 

152 # Initialize overpass connection 

153 api = overpy.Overpass() 

154 

155 # Request data from API 

156 logger.info(msg='querying OSM database for infra types...') 

157 # query_cmd = "way({},{},{},{}) [""highway=motorway_link""]; (._;>;); out body;" 

158 

159 query_cmd = "[bbox: {},{},{},{}];(" 

160 

161 if bridge_highway: 

162 logger.info(msg='\t - bridge_highway') 

163 query_cmd += 'way[highway~"^(motorway|motorway_link|trunk|trunk_link)$"][bridge];' 

164 

165 if bridge_railway: 

166 logger.info(msg='\t - bridge_railway') 

167 query_cmd += 'way[railway=rail][bridge];' 

168 

169 if (bridge_highway is False) & (bridge_railway is False): 

170 logger.info(msg='\t - all bridges') 

171 query_cmd += 'way[bridge];' 

172 

173 query_cmd += ");(._; >;); out body;" # skel 

174 

175 cmd = query_cmd.format(ll_corner_wgs[0], ll_corner_wgs[1], 

176 ur_corner_wgs[0], ur_corner_wgs[1]) 

177 logger.info(msg="\n" + cmd + "\n") 

178 result = api.query(cmd) 

179 

180 if len(result.ways) == 0: 

181 logger.error(msg='Empty OSM query results. No bridges found.') 

182 raise ValueError 

183 

184 logger.info(msg='...done.') 

185 return result