kopia lustrzana https://github.com/OpenDroneMap/ODM
122 wiersze
3.8 KiB
Python
122 wiersze
3.8 KiB
Python
import fiona
|
|
import fiona.crs
|
|
import os
|
|
import io
|
|
import json
|
|
from opendm import system
|
|
from pyproj import CRS
|
|
from opendm.location import transformer
|
|
from opendm.utils import double_quote
|
|
from osgeo import ogr
|
|
from opendm.shots import get_origin
|
|
|
|
def compute_boundary_from_shots(reconstruction_json, buffer=0, reconstruction_offset=(0, 0)):
|
|
if not os.path.isfile(reconstruction_json):
|
|
raise IOError(reconstruction_json + " does not exist.")
|
|
|
|
with open(reconstruction_json) as f:
|
|
data = json.load(f)
|
|
reconstruction = data[0]
|
|
|
|
mp = ogr.Geometry(ogr.wkbMultiPoint)
|
|
|
|
for shot_image in reconstruction['shots']:
|
|
shot = reconstruction['shots'][shot_image]
|
|
if shot.get('gps_dop', 999999) < 999999:
|
|
camera = reconstruction['cameras'][shot['camera']]
|
|
|
|
p = ogr.Geometry(ogr.wkbPoint)
|
|
origin = get_origin(shot)
|
|
|
|
p.AddPoint_2D(origin[0] + reconstruction_offset[0], origin[1] + reconstruction_offset[1])
|
|
mp.AddGeometry(p)
|
|
|
|
if mp.GetGeometryCount() < 3:
|
|
return None
|
|
|
|
convexhull = mp.ConvexHull()
|
|
boundary = convexhull.Buffer(buffer)
|
|
|
|
return load_boundary(boundary.ExportToJson())
|
|
|
|
def load_boundary(boundary_json, reproject_to_proj4=None):
|
|
if not isinstance(boundary_json, str):
|
|
boundary_json = json.dumps(boundary_json)
|
|
|
|
with fiona.open(io.BytesIO(boundary_json.encode('utf-8')), 'r') as src:
|
|
if len(src) != 1:
|
|
raise IOError("Boundary must have a single polygon (found: %s)" % len(src))
|
|
|
|
geom = src[0]['geometry']
|
|
|
|
if geom['type'] != 'Polygon':
|
|
raise IOError("Boundary must have a polygon feature (found: %s)" % geom['type'])
|
|
|
|
rings = geom['coordinates']
|
|
|
|
if len(rings) == 0:
|
|
raise IOError("Boundary geometry has no rings")
|
|
|
|
coords = rings[0]
|
|
if len(coords) == 0:
|
|
raise IOError("Boundary geometry has no coordinates")
|
|
|
|
dimensions = len(coords[0])
|
|
|
|
if reproject_to_proj4 is not None:
|
|
t = transformer(CRS.from_proj4(fiona.crs.to_string(src.crs)),
|
|
CRS.from_proj4(reproject_to_proj4))
|
|
coords = [t.TransformPoint(*c)[:dimensions] for c in coords]
|
|
|
|
return coords
|
|
|
|
def boundary_offset(boundary, reconstruction_offset):
|
|
if boundary is None or reconstruction_offset is None:
|
|
return boundary
|
|
|
|
res = []
|
|
dims = len(boundary[0])
|
|
for c in boundary:
|
|
if dims == 2:
|
|
res.append((c[0] - reconstruction_offset[0], c[1] - reconstruction_offset[1]))
|
|
else:
|
|
res.append((c[0] - reconstruction_offset[0], c[1] - reconstruction_offset[1], c[2]))
|
|
|
|
return res
|
|
|
|
def as_polygon(boundary):
|
|
if boundary is None:
|
|
return None
|
|
|
|
return "POLYGON((" + ", ".join([" ".join(map(str, c)) for c in boundary]) + "))"
|
|
|
|
def as_geojson(boundary):
|
|
return '{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[%s]}}]}' % str(list(map(list, boundary)))
|
|
|
|
def export_to_bounds_files(boundary, proj4, bounds_json_file, bounds_gpkg_file):
|
|
with open(bounds_json_file, "w") as f:
|
|
f.write(json.dumps({
|
|
"type": "FeatureCollection",
|
|
"name": "bounds",
|
|
"features": [{
|
|
"type": "Feature",
|
|
"properties": {},
|
|
"geometry": {
|
|
"type": "Polygon",
|
|
"coordinates": [boundary]
|
|
}
|
|
}]
|
|
}))
|
|
|
|
if os.path.isfile(bounds_gpkg_file):
|
|
os.remove(bounds_gpkg_file)
|
|
|
|
kwargs = {
|
|
'proj4': proj4,
|
|
'input': double_quote(bounds_json_file),
|
|
'output': double_quote(bounds_gpkg_file)
|
|
}
|
|
|
|
system.run('ogr2ogr -overwrite -f GPKG -a_srs "{proj4}" {output} {input}'.format(**kwargs))
|
|
|