Coverage for sarvey/sarvey

1#!/usr/bin/env python

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

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/>.

30"""Download openstreetmap data for area of interest."""

31import argparse

32import logging

33import os

34import sys

35import time

36from os.path import join

37import geopandas as gpd

38from shapely import ops, Point

39import matplotlib

41from sarvey.osm_utils import runOsmQueryBridge, runOsmQuery, getSpatialExtend

43try:

44 matplotlib.use('TkAgg')

45except ImportError as e:

46 print(e)

49EXAMPLE = """Example:

50 sarvey_osm --geom ./geometryRadar.h5 --railway # download railway

51 sarvey_osm --geom ./geometryRadar.h5 --highway # download highway

52 sarvey_osm --geom ./geometryRadar.h5 --railway --bridge # download railway bridge

53 sarvey_osm --geom ./geometryRadar.h5 --railway -o mask_railway.shp # specify output path

54"""

57def create_parser():

58 """Create_parser."""

59 parser = argparse.ArgumentParser(

60 description='Download transport infrastructure information from openstreetmap and store as shp-file.',

61 formatter_class=argparse.RawTextHelpFormatter,

62 epilog=EXAMPLE)

64 parser.add_argument('-w', '--work_dir', dest='work_dir', default=None,

65 help='absolute path to working directory\n' +

66 '(default: current directory).')

68 parser.add_argument('--geom', dest='geom_file', default=None,

69 help='path to existing geometryRadar.h5 file')

71 parser.add_argument('--railway', dest='railway', action="store_true", default=False,

72 help='Set true to query railways.')

74 parser.add_argument('--highway', dest='highway', action="store_true", default=False,

75 help='Set true to query highways.')

77 parser.add_argument('--bridge', dest='bridge', action="store_true", default=False,

78 help='Set true to mask bridges.\n' +

79 'If --railway or --highway set true, only railway/highway bridges are queried.')

81 parser.add_argument('-o', dest='out_file_name', default='osm_infra.shp',

82 help="name of output file. (default: 'osm_infra.shp')")

84 return parser

87def downloadOSM(*, railway: bool, highway: bool, bridge: bool,

88 work_dir: str, out_file_name: str, logger: logging.Logger, geom_file: str):

89 """Download openstreetmap data and store to file.

91 Parameters

92 ----------

93 railway: bool

94 download railway data.

95 highway: bool

96 download highway data.

97 bridge: bool

98 download bridge data.

99 work_dir: str

100 working directory.

101 out_file_name: str

102 output file name.

103 logger: logging.Logger

104 logger.

105 geom_file: str

106 path to geometryRadar.h5 file.

107 """

108 logger.info(msg="Start creating mask file based on openstreetmap data.")

109

110 # get bounding box

111 ll_bbox, ur_bbox, img_ext, coord, atr = getSpatialExtend(geom_file=geom_file, logger=logger)

112

113 # store image extend

114 gdf = gpd.GeoDataFrame({"geometry": gpd.geoseries.GeoSeries(img_ext)})

115 gdf = gdf.dissolve().convex_hull

116 gdf.to_file(join(work_dir, "img_extend.gpkg"))

117

118 # store bounding box

119 bbox_points = [

120 Point(ll_bbox[1], ll_bbox[0]),

121 Point(ur_bbox[1], ll_bbox[0]),

122 Point(ur_bbox[1], ur_bbox[0]),

123 Point(ll_bbox[1], ur_bbox[0])

124 ]

125

126 gdf = gpd.GeoDataFrame({"geometry": gpd.geoseries.GeoSeries(bbox_points)})

127 gdf = gdf.dissolve().convex_hull

128 gdf.to_file(join(work_dir, "bounding_box.gpkg"))

129

130 if (not railway) & (not highway) & (not bridge):

131 logger.error(msg="No infrastructure type was specified.")

132 return

133

134 if bridge:

135 # get requested OSM layer

136 query_result = runOsmQueryBridge(

137 ll_corner_wgs=ll_bbox, ur_corner_wgs=ur_bbox,

138 bridge_highway=highway, bridge_railway=railway,

139 logger=logger

140 )

141 else:

142 type_list = list()

143 if railway:

144 type_list += ["rail"]

145 if highway:

146 type_list += ["motorway", "motorway_link", "trunk", "trunk_link"]

147

148 # get requested OSM layer

149 query_result = runOsmQuery(ll_corner_wgs=ll_bbox, ur_corner_wgs=ur_bbox,

150 type_list=type_list, logger=logger)

151

152 multi_line_list = list()

153 for way in query_result.ways:

154 if "area" in way.tags:

155 if way.tags["area"] == "yes":

156 logger.info('Area flag is true')

157 continue

158 else:

159 # keep coordinates in lat/lon. It will be needed in masking step.

160 coord = [[float(way.nodes[i].lon), float(way.nodes[i].lat)] for i in range(len(way.nodes))]

161 multi_line_list.append(coord)

162

163 # Merge all road segments

164 merged_road = list(ops.linemerge(multi_line_list).geoms)

165 gdf = gpd.GeoDataFrame({"geometry": gpd.GeoSeries(merged_road)})

166 # gdf = gdf.set_crs(crs=utm_crs) # set appropriate CRS

167 # todo: add attributes if required

168

169 # todo: check ending of output file name

170 gdf.to_file(join(work_dir, out_file_name))

171 logger.info(msg="OSM download finished.")

172

173

174def main(iargs=None):

175 """Download openstreetmap data and store to file."""

176 # check input

177 parser = create_parser()

178 inps = parser.parse_args(args=iargs)

179

180 # initiate logger

181 logging_level = logging.getLevelName('DEBUG')

182

183 log_format = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')

184 logger = logging.getLogger(__name__)

185

186 current_datetime = time.strftime("%Y-%m-%d-%H-%M-%S", time.localtime())

187 log_filename = f"sarvey_osm_{current_datetime}.log"

188 if not os.path.exists(os.path.join(os.getcwd(), "logfiles")):

189 os.mkdir(os.path.join(os.getcwd(), "logfiles"))

190 file_handler = logging.FileHandler(filename=os.path.join(os.getcwd(), "logfiles", log_filename))

191 file_handler.setFormatter(log_format)

192 logger.addHandler(file_handler)

193

194 console_handler = logging.StreamHandler(sys.stdout)

195 console_handler.setFormatter(log_format)

196 logger.addHandler(console_handler)

197 logger.setLevel(logging_level)

198

199 if inps.work_dir is None:

200 work_dir = os.getcwd()

201 else:

202 work_dir = inps.work_dir

203 if not os.path.exists(path=work_dir):

204 logger.info(msg='create output folder: ' + work_dir)

205 os.mkdir(path=work_dir)

206 logger.info(msg='working directory: {}'.format(work_dir))

207

208 downloadOSM(

209 railway=inps.railway,

210 highway=inps.highway,

211 bridge=inps.bridge,

212 work_dir=work_dir,

213 out_file_name=inps.out_file_name,

214 logger=logger,

215 geom_file=inps.geom_file

216 )

217

218

219if __name__ == '__main__':

220 main()

Coverage for sarvey/sarvey_osm.py: 0%

82 statements