EPSG code support, GCP files reprojection to UTM

Former-commit-id: a247eb654b
pull/1161/head
Piero Toffanin 2019-06-21 13:14:44 -04:00
rodzic e774ca794c
commit 2875efea51
12 zmienionych plików z 154 dodań i 45 usunięć

Wyświetl plik

@ -95,17 +95,13 @@ def config():
choices=processopts,
help=('Can be one of:' + ' | '.join(processopts)))
parser.add_argument('--video',
metavar='<string>',
help='Path to the video file to process')
# parser.add_argument('--video',
# metavar='<string>',
# help='Path to the video file to process')
parser.add_argument('--slam-config',
metavar='<string>',
help='Path to config file for orb-slam')
parser.add_argument('--proj',
metavar='<PROJ4 string>',
help='Projection used to transform the model into geographic coordinates')
# parser.add_argument('--slam-config',
# metavar='<string>',
# help='Path to config file for orb-slam')
parser.add_argument('--min-num-features',
metavar='<integer>',

Wyświetl plik

@ -6,6 +6,9 @@ from pyproj import CRS
class GCPFile:
def __init__(self, gcp_path):
if not gcp_path:
raise RuntimeError("gcp_path must be specified")
self.gcp_path = gcp_path
self.entries = []
self.raw_srs = ""
@ -30,39 +33,64 @@ class GCPFile:
else:
log.ODM_WARNING("Malformed GCP line: %s" % line)
def entries_dict(self):
def iter_entries(self, allowed_filenames=None):
for entry in self.entries:
pe = self.parse_entry(entry)
if allowed_filenames is None or pe.filename in allowed_filenames:
yield pe
def parse_entry(self, entry):
if entry:
parts = entry.split()
x, y, z, px, py, filename = parts[:6]
extras = " ".join(parts[6:])
yield {
'x': x,
'y': y,
'z': z,
'px': px,
'py': py,
'filename': filename,
'extras': extras
}
def entry_dict_to_s(self, entry):
return "{x} {y} {z} {px} {py} {filename} {extras}".format(**entry).rstrip()
return GCPEntry(float(x), float(y), float(z), int(px), int(py), filename, extras)
def get_entry(self, n):
if n < self.entries_count():
return self.parse_entry(self.entries[n])
def entries_count(self):
return len(self.entries)
def exists(self):
return self.gcp_path and os.path.exists(self.gcp_path)
return bool(self.gcp_path and os.path.exists(self.gcp_path))
def wgs84_utm_zone(self):
"""
Finds the UTM zone where the first point of the GCP falls into
:return utm zone string valid for a coordinates header
"""
if len(self.entries) > 0:
point = list(self.entries_dict())[0]
if self.entries_count() > 0:
entry = self.get_entry(0)
longlat = CRS.from_epsg("4326")
lat, lon = location.transform(self.srs, longlat, point['x'], point['y'])
lon, lat = location.transform2(self.srs, longlat, entry.x, entry.y)
utm_zone, hemisphere = location.get_utm_zone_and_hemisphere_from(lon, lat)
return "WGS84 UTM %s%s" % (utm_zone, hemisphere)
def create_utm_copy(self, gcp_file_output, filenames=None):
"""
Creates a new GCP file from an existing GCP file
by optionally including only filenames and reprojecting each point to
a UTM CRS
"""
if os.path.exists(gcp_file_output):
os.remove(gcp_file_output)
output = [self.wgs84_utm_zone()]
target_srs = location.parse_srs_header(output[0])
transformer = location.transformer(self.srs, target_srs)
for entry in self.iter_entries(filenames):
entry.x, entry.y, entry.z = transformer.TransformPoint(entry.x, entry.y, entry.z)
output.append(str(entry))
with open(gcp_file_output, 'w') as f:
f.write('\n'.join(output) + '\n')
return gcp_file_output
def make_filtered_copy(self, gcp_file_output, images_dir, min_images=3):
"""
Creates a new GCP file from an existing GCP file includes
@ -81,9 +109,9 @@ class GCPFile:
output = [self.raw_srs]
files_found = 0
for entry in self.entries_dict():
if entry['filename'] in files:
output.append(self.entry_dict_to_s(entry))
for entry in self.iter_entries():
if entry.filename in files:
output.append(str(entry))
files_found += 1
if files_found >= min_images:
@ -91,3 +119,19 @@ class GCPFile:
f.write('\n'.join(output) + '\n')
return gcp_file_output
class GCPEntry:
def __init__(self, x, y, z, px, py, filename, extras=""):
self.x = x
self.y = y
self.z = z
self.px = px
self.py = py
self.filename = filename
self.extras = extras
def __str__(self):
return "{} {} {} {} {} {} {}".format(self.x, self.y, self.z,
self.px, self.py,
self.filename,
self.extras).rstrip()

Wyświetl plik

@ -1,6 +1,7 @@
import math
from opendm import log
from pyproj import Proj, Transformer, CRS
from osgeo import osr
def extract_utm_coords(photos, images_path, output_coords_file):
"""
@ -54,10 +55,32 @@ def extract_utm_coords(photos, images_path, output_coords_file):
for coord in coords:
f.write("%s %s %s\n" % (coord[0] - dx, coord[1] - dy, coord[2]))
def transform(from_srs, to_srs, x, y):
transformer = Transformer.from_crs(from_srs, to_srs)
return transformer.transform(x, y)
def transform2(from_srs, to_srs, x, y):
return transformer(from_srs, to_srs).TransformPoint(x, y, 0)[:2]
def transform3(from_srs, to_srs, x, y, z):
return transformer(from_srs, to_srs).TransformPoint(x, y, z)
def proj_srs_convert(srs):
"""
Convert a Proj SRS object to osr SRS object
"""
res = osr.SpatialReference()
epsg = srs.to_epsg()
if epsg:
res.ImportFromEPSG(epsg)
else:
proj4 = srs.to_proj4()
res.ImportFromProj4(proj4)
return res
def transformer(from_srs, to_srs):
src = proj_srs_convert(from_srs)
tgt = proj_srs_convert(to_srs)
return osr.CoordinateTransformation(src, tgt)
def get_utm_zone_and_hemisphere_from(lon, lat):
"""
Calculate the UTM zone and hemisphere that a longitude/latitude pair falls on
@ -107,7 +130,7 @@ def parse_srs_header(header):
'datum': datum
}
proj4 = '+proj=utm +zone={zone} +datum={datum} +no_defs=True'
proj4 = '+proj=utm +zone={zone} +datum={datum} +units=m +no_defs=True'
if utm_pole == 'S':
proj4 += ' +south=True'

Wyświetl plik

@ -58,7 +58,7 @@ class OSFMContext:
raise Exception("Reconstruction could not be generated")
def setup(self, args, images_path, photos, gcp_path=None, append_config = [], rerun=False):
def setup(self, args, images_path, photos, gcp_path=None, append_config = [], rerun=False): #georeferenced=True,
"""
Setup a OpenSfM project
"""
@ -116,6 +116,10 @@ class OSFMContext:
else:
config.append("align_method: orientation_prior")
config.append("align_orientation_prior: vertical")
# if not georeferenced:
# config.append("")
if args.use_hybrid_bundle_adjustment:
log.ODM_DEBUG("Enabling hybrid bundle adjustment")

Wyświetl plik

@ -103,7 +103,7 @@ class ODMLoadDatasetStage(types.ODM_Stage):
reconstruction = types.ODM_Reconstruction(photos)
if tree.odm_georeferencing_gcp:
reconstruction.georeference_with_gcp(tree.odm_georeferencing_gcp, tree.odm_georeferencing_coords)
reconstruction.georeference_with_gcp(tree.odm_georeferencing_gcp, tree.odm_georeferencing_coords, reload_coords=self.rerun())
else:
reconstruction.georeference_with_gps(tree.dataset_raw, tree.odm_georeferencing_coords, reload_coords=self.rerun())

1
tests/assets/.gitignore vendored 100644
Wyświetl plik

@ -0,0 +1 @@
output/

Wyświetl plik

@ -1,3 +0,0 @@
EPSG:4326
44.701151842605441 -85.61322729531112 150.5 1331 350 DJI_0002.JPG
44.701023803183617 -85.612848642726888 151.5 1331 350 DJI_0003.JPG

Wyświetl plik

@ -0,0 +1,3 @@
EPSG:4326
105.03506 -5.23270 20 1331 350 DJI_0002.JPG

Wyświetl plik

@ -0,0 +1,3 @@
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
-85.61322729531112 44.701151842605441 150.5 1331 350 DJI_0002.JPG
-85.612848642726888 44.701023803183617 151.5 1331 350 DJI_0003.JPG

Wyświetl plik

@ -0,0 +1,2 @@
EPSG:2251
26607594.75468351 -26813.882366498266 563.1987598425 1331 350 DJI_0002.JPG

Wyświetl plik

@ -1,21 +1,57 @@
import time
import unittest
import os
import shutil
from opendm.gcp import GCPFile
class TestRemote(unittest.TestCase):
class TestGcp(unittest.TestCase):
def setUp(self):
pass
if os.path.exists("tests/assets/output"):
shutil.rmtree("tests/assets/output")
os.makedirs("tests/assets/output")
# pass
def test_utm_north(self):
gcp = GCPFile("tests/assets/gcp_utm_north_valid.txt")
self.assertTrue(gcp.exists())
self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 16N")
def test_epsg(self):
gcp = GCPFile("tests/assets/gcp_epsg_valid.txt")
def test_latlon_south(self):
gcp = GCPFile("tests/assets/gcp_latlon_south.txt")
self.assertTrue(gcp.exists())
self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 48S")
def test_latlon(self):
gcp = GCPFile("tests/assets/gcp_latlon_valid.txt")
self.assertTrue(gcp.exists())
self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 16N")
def test_utm_conversion(self):
gcp = GCPFile("tests/assets/gcp_latlon_valid.txt")
copy = GCPFile(gcp.create_utm_copy("tests/assets/output/gcp_utm.txt"))
self.assertTrue(copy.exists())
self.assertEqual(copy.raw_srs, "WGS84 UTM 16N")
self.assertEqual(copy.get_entry(0).x, 609865.707705)
self.assertEqual(copy.get_entry(0).y, 4950688.36182)
def test_utm_conversion_feet(self):
gcp = GCPFile("tests/assets/gcp_michigan_feet_valid.txt")
copy = GCPFile(gcp.create_utm_copy("tests/assets/output/gcp_utm_z.txt"))
self.assertTrue(copy.exists())
self.assertEqual(copy.raw_srs, "WGS84 UTM 16N")
self.assertEqual(round(copy.get_entry(0).x, 3), 609925.818)
self.assertEqual(round(copy.get_entry(0).y, 3), 4950688.772)
self.assertEqual(round(copy.get_entry(0).z, 3), 171.663)
def test_filtered_copy(self):
gcp = GCPFile('tests/assets/gcp_latlon_valid.txt')
self.assertTrue(gcp.exists())
self.assertEqual(gcp.entries_count(), 2)
copy = GCPFile(gcp.make_filtered_copy('tests/assets/output/filtered_copy.txt', 'tests/assets/images', min_images=1))
self.assertTrue(copy.exists())
self.assertEqual(copy.entries_count(), 1)
if __name__ == '__main__':
unittest.main()