Coverage for aisdb/gis.py: 98%

1'''Geometry and GIS utilities'''

3import os

4import pathlib

5import tempfile

6from datetime import datetime, timedelta

7from functools import partial

9import numpy as np

10import shapely.ops

11import shapely.geometry

12import warnings

13from shapely.geometry import Polygon, LineString, Point

15from aisdb.aisdb import haversine

16from aisdb.proc_util import glob_files

19def shiftcoord(x, rng=180):

20 ''' Correct longitude coordinates to be within range(-180, 180)

21 using a linear shift and modulus.

22 For latitude coordinate correction, set rng to 90.

24 For example: longitude 181 would be corrected to -179 deg.

25 '''

26 assert len(x) > 0, 'x must be array-like'

27 if not isinstance(x, np.ndarray):

28 x = np.array(x)

29 shift_idx = np.where(np.abs(x) != rng)[0]

30 for idx in shift_idx:

31 x[idx] = ((x[idx] + rng) % 360) - rng

32 flip_idx = np.where(np.abs(x) == rng)[0]

33 for idx in flip_idx:

34 x[idx] *= -1

35 assert (rng * -1 <= np.all(x) <= rng)

36 return x

39def dt_2_epoch(dt_arr, t0=datetime(1970, 1, 1, 0, 0, 0)):

40 ''' convert datetime.datetime to epoch minutes '''

41 delta = lambda dt: (dt - t0).total_seconds()

42 if isinstance(dt_arr, (list, np.ndarray)):

43 return np.array(list(map(float, map(delta, dt_arr))))

44 elif isinstance(dt_arr, (datetime)):

45 return int(delta(dt_arr))

46 else:

47 raise ValueError('input must be datetime or array of datetimes')

50def epoch_2_dt(ep_arr, t0=datetime(1970, 1, 1, 0, 0, 0), unit='seconds'):

51 ''' convert epoch minutes to datetime.datetime '''

53 delta = lambda ep, unit: t0 + timedelta(**{unit: ep})

55 if isinstance(ep_arr, (list, np.ndarray)):

56 return np.array(list(map(partial(delta, unit=unit), map(int, ep_arr))))

58 elif isinstance(ep_arr,

59 (float, int, np.uint32, np.int32, np.uint64, np.int64)):

60 return delta(int(ep_arr), unit=unit)

62 else:

63 raise ValueError(

64 f'input must be integer or array of integers. got {ep_arr=}{type(ep_arr)}'

65 )

68def delta_meters(track, rng=None):

69 ''' compute haversine distance in meters between track positions for a

70 given track

72 args:

73 track (dict)

74 track vector dictionary

75 rng (range)

76 optionally restrict computed values to given index range

77 '''

78 rng = range(len(track['time'])) if rng is None else rng

79 return np.array(

80 list(

81 map(haversine, track['lon'][rng][:-1], track['lat'][rng][:-1],

82 track['lon'][rng][1:], track['lat'][rng][1:])))

85def delta_seconds(track, rng=None):

86 ''' compute elapsed time between track positions for a given track

88 args:

89 track (dict)

90 track vector dictionary

91 rng (range)

92 optionally restrict computed values to given index range

93 '''

94 if isinstance(track['time'], list):

95 track['time'] = np.array(track['time'])

96 assert isinstance(track['time'], np.ndarray), f'got {track["time"] = }'

97 rng = range(len(track['time'])) if rng is None else rng

98 return np.array(list((track['time'][rng][1:] - track['time'][rng][:-1])))

100

101def delta_knots(track, rng=None):

102 ''' compute speed over ground in knots between track positions for a given

103 track using (haversine distance / time)

104

105 args:

106 track (dict)

107 track vector dictionary

108 rng (range)

109 optionally restrict computed values to given index range

110 '''

111 rng = range(len(track['time'])) if rng is None else rng

112 ds = np.array([np.max((1, s)) for s in delta_seconds(track, rng)],

113 dtype=object)

114 return delta_meters(track, rng) / ds * 1.9438445

115

116

117def radial_coordinate_boundary(x, y, radius=100000):

118 ''' Defines a bounding box area for a given point and radial

119 distance in meters. Returns degree boundaries with a minimum diameter

120 of approximately 2 * ``radius`` meters.

121

122 The boundaries are approximated by converting input coordinates

123 from degrees to radians, and computing a radial delta by

124 dividing an input value by the earth radius. The radial delta

125 is added or subtracted from the input point for each cardinal

126 direction, and then converted back from radians to degrees.

127

128 args:

129 x (float)

130 longitude

131 y (float)

132 latitude

133 radius (int, float)

134 minimum radial distance

135 '''

136 # radians

137 earth_radius_m = 6371088

138 rlon = np.pi * x / 180

139 rlat = np.pi * y / 180

140 parallel_radius = earth_radius_m * np.cos(rlat)

141

142 # radial delta

143 rlonmin = rlon - radius / parallel_radius

144 rlonmax = rlon + radius / parallel_radius

145 rlatmin = rlat - radius / earth_radius_m

146 rlatmax = rlat + radius / earth_radius_m

147

148 return {

149 'xmin': 180 * rlonmin / np.pi,

150 'xmax': 180 * rlonmax / np.pi,

151 'ymin': 180 * rlatmin / np.pi,

152 'ymax': 180 * rlatmax / np.pi

153 }

154

155

156def distance3D(x1, y1, x2, y2, depth_metres):

157 ''' haversine/pythagoras approximation of vessel distance to

158 point at given depth

159 '''

160 a2 = haversine(x1=x1, y1=y1, x2=x2, y2=y2)**2

161 b2 = abs(depth_metres)**2

162 c2 = a2 + b2

163 return np.sqrt(c2)

164

165

166def vesseltrack_3D_dist(tracks, x1, y1, z1, colname='distance_metres'):

167 ''' appends approximate distance to a submerged point at every

168 surface-level position. distance is approximated using the haversine

169 function and pythagoras.

170

171 x1 (float)

172 point longitude

173 y1 (float)

174 point latitude

175 z1 (float)

176 point depth (metres)

177 colname (string)

178 track dictionary key for which depth values will be set.

179 by default, distances are appended to the 'distance_metres'

180 key

181 '''

182 for track in tracks:

183 track['dynamic'] = track['dynamic'].union(set([colname]))

184 dists = [

185 distance3D(x1=x1, y1=y1, x2=x, y2=y, depth_metres=z1)

186 for x, y in zip(track['lon'], track['lat'])

187 ]

188 track[colname] = np.array(dists, dtype=object)

189 yield track

190

191

192def mask_in_radius_2D(tracks, xy, distance_meters):

193 ''' radial filtering using great circle distance at surface level

194

195 tracks (:class:`aisdb.track_gen.TrackGen`)

196 track dictionary iterator

197 xy (tuple of floats)

198 target longitude and latitude coordinate pair

199 distance_meters (int)

200 maximum distance in meters

201 '''

202 for track in tracks:

203 mask = [

204 haversine(track['lon'][i], track['lat'][i], xy[0], xy[1]) <

205 distance_meters for i in range(len(track['time']))

206 ]

207 if sum(mask) == 0:

208 continue

209 yield dict(

210 **{k: track[k]

211 for k in track['static']},

212 **{k: track[k][mask]

213 for k in track['dynamic']},

214 static=track['static'],

215 dynamic=track['dynamic'],

216 )

217

218

219class Domain():

220 ''' collection of zone geometry dictionaries, with additional

221 statistics such as hull bounding box

222

223 args:

224 name: string

225 Domain name

226 zones: list of dictionaries

227 Collection of zone geometry dictionaries.

228 Must have keys 'name' (string) and 'geometry'

229 (shapely.geometry.Polygon)

230

231 >>> domain = Domain(name='example', zones=[{

232 ... 'name': 'zone1',

233 ... 'geometry': shapely.geometry.Polygon([(-40,60), (-40, 61), (-41, 61), (-41, 60), (-40, 60)])

234 ... }, ])

235

236 attr:

237 self.name

238 self.zones

239 self.boundary

240 self.minX

241 self.minY

242 self.maxX

243 self.maxY

244

245 '''

246

247 _meridian = LineString(

248 np.array((

249 (-180, -180, 180, 180),

250 (-90, 90, 90, -90),

251 )).T)

252

253 def _zone_max_radius(self, geom, zone_x, zone_y):

254 ''' computes the maximum distance to the centroid '''

255 return np.max([

256 haversine(geom.centroid.x, geom.centroid.y, x2, y2)

257 for x2, y2 in zip(zone_x, zone_y)

258 ])

259

260 def _add_zone(self, name, x, y):

261 '''

262 if name[-2:] == '_b':

263 if (x0b := np.min(x)) < self.minX_b:

264 self.minX_b = x0b

265 if (x0c := np.min(x)) > self.minX:

266 self.minX = x0c

267

268 elif name[-2:] == '_c':

269 if (x1b := np.max(x)) < self.maxX_c:

270 self.maxX_c = x1b

271 if (x1c := np.max(x)) > self.maxX:

272 self.maxX = x1c

273

274 else:

275 '''

276 if ((x0a := np.min(x)) < self.minX):

277 self.minX = x0a

278 if ((x1a := np.max(x)) > self.maxX):

279 self.maxX = x1a

280

281 if np.min(y) < self.minY:

282 self.minY = np.min(y)

283 if np.max(y) > self.maxY:

284 self.maxY = np.max(y)

285 '''

286 if np.min(x) < self.minX:

287 self.minX = np.min(x)

288 if np.max(x) > self.maxX:

289 self.maxX = np.max(x)

290 '''

291

292 assert self.minX < self.maxX

293 assert self.minY < self.maxY

294 assert -180 <= self.minX <= 180 and -180 <= self.maxX <= 180

295 assert -90 <= self.minY <= 90 and -90 <= self.maxY <= 90

296

297 geom = Polygon(zip(x, y))

298

299 self.zones.update({

300 name: {

301 'geometry': geom,

302 'maxradius': self._zone_max_radius(geom, x, y)

303 }

304 })

305

306 return

307

308 def _handle_outofbounds_zone(self, zone, zones_dir):

309 zones_west = zones_dir / 'west'

310 zones_east = zones_dir / 'east'

311 zones_corr = zones_dir / 'corr'

312 stringify = lambda x, y: map(

313 ','.join, zip(map(str, x), map(lambda y: y + '\n', map(str, y))))

314

315 for g in self.split_geom(zone):

316 if g.centroid.x < -180:

317 x, y = np.array(g.boundary.coords.xy)

318 if not os.path.isdir(zones_west): 318 ↛ 320line 318 didn't jump to line 320, because the condition on line 318 was never false

319 os.mkdir(zones_west)

320 with open(os.path.join(zones_west, zone['name'] + '_west.txt'),

321 'w') as w:

322 w.writelines(stringify(shiftcoord(x), y))

323 elif g.centroid.x > 180:

324 x, y = np.array(g.boundary.coords.xy)

325 if not os.path.isdir(zones_east): 325 ↛ 327line 325 didn't jump to line 327, because the condition on line 325 was never false

326 os.mkdir(zones_east)

327 with open(os.path.join(zones_east, zone['name'] + '_east.txt'),

328 'w') as w:

329 w.writelines(stringify(shiftcoord(x), y))

330 else:

331 x, y = np.array(g.boundary.coords.xy)

332 if not os.path.isdir(zones_corr):

333 os.mkdir(zones_corr)

334 with open(os.path.join(zones_corr, zone['name'] + '_corr.txt'),

335 'w') as w:

336 w.writelines(stringify(x, y))

337

338 def __init__(self, name, zones=[], **kw):

339 ''' Initialize the domain from zone geometries '''

340

341 if len(zones) == 0:

342 raise ValueError(

343 'domain needs to have atleast one polygon geometry')

344 self.name = name

345 self.zones = {}

346 self.minX, self.maxX = 180, -180

347 # self.minX_b = 180

348 self.minY, self.maxY = 90, -90

349 # self.maxX_c = -180

350

351 valid_domain = True

352 zones_dir = pathlib.Path(tempfile.mkdtemp(prefix='aisdb_zones_'))

353

354 for zone in zones:

355 if 'name' not in zone.keys():

356 raise KeyError(f'Zone missing \'name\' key: {zone=}')

357 if 'geometry' not in zone.keys():

358 raise KeyError(f'Zone missing \'geometry\' key: {zone=}')

359

360 x, y = zone['geometry'].boundary.coords.xy

361 if not (np.min(x) >= -180 and np.max(x) <= 180):

362 #warnings.warn(f'dividing geometry... {zone["name"]}')

363 valid_domain = False

364 self._handle_outofbounds_zone(zone, zones_dir)

365 else:

366 self._add_zone(zone['name'], x, y)

367

368 if not valid_domain:

369 '''

370 if not os.path.isdir(os.path.dirname(

371 zones_corr)) or not os.path.isdir(zones_corr):

372 print('Creating new output directory in ',

373 os.path.dirname(zones_dir))

374 os.makedirs(zones_dir, exist_ok=True)

375 '''

376

377 for zonename, zone in self.zones.items():

378 zone['name'] = zonename

379 self._handle_outofbounds_zone(zone, zones_dir)

380

381 raise ValueError(

382 'Invalid zone geometry! '

383 'Exceeds longitude range -180 to 180. '

384 'If you want to query a bounding box spanning 180 degrees '

385 'longitude, consider querying multiple times instead.\n'

386 f'Saved modified geometries in {str(zones_dir)}, try using these corrected domains:\n'

387 f'\tdomain1 = aisdb.DomainFromTxts(domainName=\'{name}_corr\', folder={str(zones_dir)}{os.path.sep}corrected)\n'

388 f'\tdomain2 = aisdb.DomainFromTxts(domainName=\'{name}_west\', folder={str(zones_dir)}{os.path.sep}west))\n'

389 f'\tdomain3 = aisdb.DomainFromTxts(domainName=\'{name}_east\', folder={str(zones_dir)}{os.path.sep}east))\n'

390 )

391

392 assert self.minX < self.maxX

393 assert self.minY < self.maxY

394

395 self.boundary = {

396 'xmin': self.minX,

397 'xmax': self.maxX,

398 'ymin': self.minY,

399 'ymax': self.maxY

400 }

401

402 def nearest_polygons_to_point(self, x, y):

403 ''' compute great circle distance for this point to each polygon

404 centroid, subtracting the maximum polygon radius.

405 returns all zones with distances less than zero meters, sorted by

406 nearest first

407 '''

408 assert float(x) or x == 0.0, f'{type(x)} {x=}{y=}'

409 assert float(y) or y == 0.0, f'{type(y)} {x=}{y=}'

410 assert isinstance(self.zones, dict)

411 dist_to_centroids = {}

412 for name, z in self.zones.items():

413 dist_to_centroids.update({

414 name:

415 haversine(

416 x,

417 y,

418 z['geometry'].centroid.x,

419 z['geometry'].centroid.y,

420 ) - z['maxradius']

421 })

422 return dist_to_centroids

423

424 def point_in_polygon(self, x, y):

425 ''' Returns the zone containing the given coordinates.

426 if there are multiple zones containing the coordinates,

427 the zone with the nearest centroid will be selected.

428

429 args:

430 x (float)

431 longitude value

432 y (float)

433 latitude value

434 '''

435 assert float(x) or x == 0.0, f'{type(x)} {x=}{y=}'

436 assert float(y) or y == 0.0, f'{type(y)} {x=}{y=}'

437 assert len(self.zones) > 0

438 # first pass filter using distance to centroid, subtracting max radius.

439 # discard all geometry with a distance over zero

440 nearest = {

441 k: v

442 for k, v in sorted(self.nearest_polygons_to_point(x, y).items(),

443 key=lambda item: item[1]) if v < 0

444 }

445 # check for zone containment, starting with the nearest centroid

446 for key, value in nearest.items():

447 if self.zones[key]['geometry'].contains(Point(x, y)):

448 return key

449 return 'Z0'

450

451 def split_geom(self, zone):

452 ''' Ensure that the zone doesn't intersect longitude 180 or -180.

453 If it does, divide it into two zones.

454 '''

455 merged = shapely.ops.linemerge(

456 [zone['geometry'].boundary, self._meridian])

457 border = shapely.ops.unary_union(merged)

458 decomp = shapely.ops.polygonize(border)

459 return decomp

460

461

462class DomainFromTxts(Domain):

463 ''' subclass of :class:`aisdb.gis.Domain`. used for convenience to load

464 zone geometry from .txt files directly

465 '''

466

467 def __init__(self, domainName, folder, ext='txt'):

468

469 files = glob_files(folder, ext=ext)

470 zones = []

471 for txt in files:

472 filename = txt.rsplit(os.path.sep, 1)[1].split('.')[0]

473 with open(txt, 'r') as f:

474 pts = f.read()

475 xy = list(

476 map(float,

477 pts.replace('\n\n', '').replace('\n',

478 ',').split(',')[:-1]))

479 x, y = np.array(xy[::2]), np.array(xy[1::2])

480 geom = Polygon(zip(x, y))

481 zones.append({'name': filename, 'geometry': geom})

482 super().__init__(domainName, zones, setattrs=False)

483

484

485class DomainFromPoints(Domain):

486 ''' Subclass of :class:`aisdb.gis.Domain`. Used for convenience to generate

487 bounding box polygons from longitude/latitude pairs and radial

488 distances, where the minimum radius is specified in meters.

489 '''

490

491 def __init__(self,

492 points,

493 radial_distances,

494 names=[],

495 domainName='domain'):

496 ''' Creates bounding-box polygons having a minimum radius atleast

497 approximately ``radial_distances`` in metres, centered on

498 ``points``.

499

500 args:

501 points (list)

502 coordinate XY pairs

503 radial_distances (list)

504 approximate distance in meters to extend the bounding box.

505 the distance given will be used as the minimum distance to

506 box boundaries

507 names (list)

508 optionally assign a zone name for each point

509

510

511 '''

512 if names == []: 512 ↛ 514line 512 didn't jump to line 514, because the condition on line 512 was never false

513 names = [f'{i:02d}' for i in range(len(points))]

514 zones = []

515 for xy, d, i in zip(points, radial_distances, names):

516 bounds = radial_coordinate_boundary(*xy, d)

517 geom = {

518 'name':

519 i,

520 'geometry':

521 Polygon([

522 (bounds['xmin'], bounds['ymin']),

523 (bounds['xmin'], bounds['ymax']),

524 (bounds['xmax'], bounds['ymax']),

525 (bounds['xmax'], bounds['ymin']),

526 (bounds['xmin'], bounds['ymin']),

527 ]),

528 }

529 zones.append(geom)

530 super().__init__(name=domainName, zones=zones)