mirror
/
OpenDroneMap-ODM
kopia lustrzana https://github.com/OpenDroneMap/ODM
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

Merge master

pull/1283/head
Piero Toffanin 2021-05-11 10:49:52 -04:00
rodzic 6acc9d4162 f0c495eedc
commit 717b6dcb6e
8 zmienionych plików z 145 dodań i 860 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

2
VERSION
Wyświetl plik

@ -1 +1 @@
2.5.0
2.5.0

177
opendm/cutline.py
Wyświetl plik

@ -1,27 +1,46 @@
import os
import shutil
import rasterio
import fiona
import numpy as np
import math
from opendm import log
from opendm import io
from opendm import concurrency
from opendm import get_image_size
from opendm import system
import math
def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrency=1, tmpdir=None, scale=1):
from skimage.feature import canny
from skimage.draw import line
from skimage.graph import route_through_array
from shapely.geometry import LineString, mapping, shape
from shapely.ops import polygonize, unary_union
def write_raster(data, file):
profile = {
'driver': 'GTiff',
'width': data.shape[1],
'height': data.shape[0],
'count': 1,
'dtype': 'float32',
'transform': None,
'nodata': None,
'crs': None
}
with rasterio.open(file, 'w', **profile) as wout:
wout.write(data, 1)
def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrency=1, scale=1):
if io.file_exists(orthophoto_file) and io.file_exists(crop_area_file):
from opendm.grass_engine import grass
log.ODM_INFO("Computing cutline")
if tmpdir and not io.dir_exists(tmpdir):
system.mkdir_p(tmpdir)
scale = max(0.0001, min(1, scale))
scaled_orthophoto = None
if scale < 1:
log.ODM_INFO("Scaling orthophoto to %s%% to compute cutline" % (scale * 100))
scaled_orthophoto = os.path.join(tmpdir, os.path.basename(io.related_file_path(orthophoto_file, postfix=".scaled")))
scaled_orthophoto = io.related_file_path(orthophoto_file, postfix=".scaled")
# Scale orthophoto before computing cutline
system.run("gdal_translate -outsize {}% 0 "
"-co NUM_THREADS={} "
@ -33,36 +52,126 @@ def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrenc
orthophoto_file,
scaled_orthophoto
))
orthophoto_file = scaled_orthophoto
try:
ortho_width,ortho_height = get_image_size.get_image_size(orthophoto_file, fallback_on_error=False)
log.ODM_INFO("Orthophoto dimensions are %sx%s" % (ortho_width, ortho_height))
number_lines = int(max(8, math.ceil(min(ortho_width, ortho_height) / 256.0)))
except:
log.ODM_INFO("Cannot compute orthophoto dimensions, setting arbitrary number of lines.")
number_lines = 32
log.ODM_INFO("Number of lines: %s" % number_lines)
# open raster
f = rasterio.open(orthophoto_file)
rast = f.read(1) # First band only
height, width = rast.shape
number_lines = int(max(8, math.ceil(min(width, height) / 256.0)))
line_hor_offset = int(width / number_lines)
line_ver_offset = int(height / number_lines)
gctx = grass.create_context({'auto_cleanup' : False, 'tmpdir': tmpdir})
gctx.add_param('orthophoto_file', orthophoto_file)
gctx.add_param('crop_area_file', crop_area_file)
gctx.add_param('number_lines', number_lines)
gctx.add_param('max_concurrency', max_concurrency)
gctx.add_param('memory', int(concurrency.get_max_memory_mb(300)))
gctx.set_location(orthophoto_file)
if line_hor_offset <= 2 or line_ver_offset <= 2:
log.ODM_WARNING("Cannot compute cutline, orthophoto is too small (%sx%spx)" % (width, height))
return
cutline_file = gctx.execute(os.path.join("opendm", "grass", "compute_cutline.grass"))
if cutline_file != 'error':
if io.file_exists(cutline_file):
shutil.move(cutline_file, destination)
log.ODM_INFO("Generated cutline file: %s --> %s" % (cutline_file, destination))
gctx.cleanup()
return destination
crop_f = fiona.open(crop_area_file, 'r')
if len(crop_f) == 0:
log.ODM_WARNING("Crop area is empty, cannot compute cutline")
return
crop_poly = shape(crop_f[1]['geometry'])
crop_f.close()
linestrings = []
# Compute canny edges on first band
edges = canny(rast)
def compute_linestrings(direction):
log.ODM_INFO("Computing %s cutlines" % direction)
# Initialize cost map
cost_map = np.full((height, width), 1, dtype=np.float32)
# Write edges to cost map
cost_map[edges==True] = 0 # Low cost
# Write "barrier, floor is lava" costs
if direction == 'vertical':
lines = [((i, 0), (i, height - 1)) for i in range(line_hor_offset, width - line_hor_offset, line_hor_offset)]
points = []
pad_x = int(line_hor_offset / 2.0)
for i in range(0, len(lines)):
a,b = lines[i]
points.append(((a[0] - pad_x , a[1]), (b[0] - pad_x, b[1])))
a,b = lines[-1]
points.append(((a[0] + pad_x , a[1]), (b[0] + pad_x, b[1])))
else:
log.ODM_WARNING("Unexpected script result: %s. No cutline file has been generated." % cutline_file)
else:
log.ODM_WARNING("Could not generate orthophoto cutline. An error occured when running GRASS. No orthophoto will be generated.")
lines = [((0, j), (width - 1, j)) for j in range(line_ver_offset, height - line_ver_offset, line_ver_offset)]
points = []
pad_y = int(line_ver_offset / 2.0)
for i in range(0, len(lines)):
a,b = lines[i]
points.append(((a[0] , a[1] - pad_y), (b[0], b[1] - pad_y)))
a,b = lines[-1]
points.append(((a[0] , a[1] + pad_y), (b[0], b[1] + pad_y)))
for a, b in lines:
rr,cc = line(*a, *b)
cost_map[cc, rr] = 9999 # Lava
# Calculate route
for a, b in points:
line_coords, cost = route_through_array(cost_map, (a[1], a[0]), (b[1], b[0]), fully_connected=True, geometric=True)
# Convert to geographic
geo_line_coords = [f.xy(*c) for c in line_coords]
# Simplify
ls = LineString(geo_line_coords)
linestrings.append(ls.simplify(0.05, preserve_topology=False))
compute_linestrings('vertical')
compute_linestrings('horizontal')
# Generate polygons and keep only those inside the crop area
log.ODM_INFO("Generating polygons... this could take a bit.")
polygons = []
for p in polygonize(unary_union(linestrings)):
if crop_poly.contains(p):
polygons.append(p)
# This should never happen
if len(polygons) == 0:
log.ODM_WARNING("No polygons, cannot compute cutline")
return
log.ODM_INFO("Merging polygons")
cutline_polygons = unary_union(polygons)
largest_cutline = cutline_polygons[0]
max_area = largest_cutline.area
for p in cutline_polygons:
if p.area > max_area:
max_area = p.area
largest_cutline = p
log.ODM_INFO("Largest cutline found: %s m^2" % max_area)
meta = {
'crs': {'init': str(f.crs).lower() },
'driver': 'GPKG',
'schema': {
'properties': {},
'geometry': 'Polygon'
}
}
# Remove previous
if os.path.exists(destination):
os.remove(destination)
with fiona.open(destination, 'w', **meta) as sink:
sink.write({
'geometry': mapping(largest_cutline),
'properties': {}
})
f.close()
log.ODM_INFO("Wrote %s" % destination)
# Cleanup
if scaled_orthophoto is not None and os.path.exists(scaled_orthophoto):
os.remove(scaled_orthophoto)
else:
log.ODM_WARNING("We've been asked to compute cutline, but either %s or %s is missing. Skipping..." % (orthophoto_file, crop_area_file))

641
opendm/grass/addons/i.cutlinesmod.py
Wyświetl plik

@ -1,641 +0,0 @@
#!/usr/bin/env python3
############################################################################
#
# MODULE: i.cutlinesmod
# AUTHOR(S): Moritz Lennert, with help of Stefanos Georganos, modified by
# Piero Toffanin
#
# PURPOSE: Create tiles the borders of which do not cut across semantically
# meaningful objects
# COPYRIGHT: (C) 1997-2018 by the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#############################################################################
#%Module
#% description: Creates semantically meaningful tile borders
#% keyword: imagery
#% keyword: tiling
#%end
#
#%option G_OPT_R_INPUT
#% description: Raster map to use as input for tiling
#% required: yes
#%end
#
#%option G_OPT_V_OUTPUT
#% description: Name of output vector map with cutline polygons
#%end
#
#%option
#% key: number_lines
#% type: integer
#% description: Number of tile border lines in each direction
#% required: yes
#%end
#
#%option
#% key: edge_detection
#% type: string
#% description: Edge detection algorithm to use
#% options: zc,canny
#% answer: zc
#% required: yes
#%end
#
#%option G_OPT_V_INPUTS
#% key: existing_cutlines
#% label: Input vector maps with existing cutlines
#% required: no
#%end
#
#%option
#% key: no_edge_friction
#% type: integer
#% description: Additional friction for non-edge pixels
#% required: yes
#% answer: 5
#%end
#
#%option
#% key: lane_border_multiplier
#% type: integer
#% description: Multiplier for borders of lanes compared to non-edge pixels
#% required: yes
#% answer: 10
#%end
#
#%option
#% key: min_tile_size
#% type: integer
#% description: Minimum size of tiles in map units
#% required: no
#%end
#
#%option
#% key: zc_threshold
#% type: double
#% label: Sensitivity of Gaussian filter (i.zc)
#% answer: 1
#% required: no
#% guisection: Zero-crossing
#%end
#
#%option
#% key: zc_width
#% type: integer
#% label: x-y extent of the Gaussian filter (i.zc)
#% answer: 9
#% required: no
#% guisection: Zero-crossing
#%end
#
#%option
#% key: canny_low_threshold
#% type: double
#% label: Low treshold for edges (i.edge)
#% answer: 3
#% required: no
#% guisection: Canny
#%end
#
#%option
#% key: canny_high_threshold
#% type: double
#% label: High treshold for edges (i.edge)
#% answer: 10
#% required: no
#% guisection: Canny
#%end
#
#%option
#% key: canny_sigma
#% type: double
#% label: Kernel radius (i.edge)
#% answer: 2
#% required: no
#% guisection: Canny
#%end
#
#%option
#% key: tile_width
#% type: integer
#% description: Width of tiles for tiled edge detection (pixels)
#% required: no
#% guisection: Parallel processing
#%end
#
#%option
#% key: tile_height
#% type: integer
#% description: Height of tiles for tiled edge detection (pixels)
#% required: no
#% guisection: Parallel processing
#%end
#
#%option
#% key: overlap
#% type: integer
#% description: Overlap between tiles for tiled edge detection (pixels)
#% required: no
#% answer: 1
#% guisection: Parallel processing
#%end
#
#%option
#% key: processes
#% type: integer
#% description: Number of parallel processes
#% answer: 1
#% required: yes
#% guisection: Parallel processing
#%end
#
#%option
#% key: memory
#% type: integer
#% description: RAM memory available (in MB)
#% answer: 300
#% required: yes
#%end
#
#%rules
#% collective: tile_width, tile_height, overlap
#%end
import os
import atexit
import grass.script as gscript
from grass.pygrass.modules.grid.grid import GridModule
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector import geometry as geom
def cleanup():
gscript.message(_("Erasing temporary files..."))
for temp_map, maptype in temp_maps:
if gscript.find_file(temp_map, element=maptype)['name']:
gscript.run_command('g.remove', flags='f', type=maptype,
name=temp_map, quiet=True)
def listzip(input1, input2):
# python3 compatible
out = zip(input1, input2)
if not isinstance(out, list):
out = list(zip(input1, input2))
return out
def main():
inputraster = options['input']
number_lines = int(options['number_lines'])
edge_detection_algorithm = options['edge_detection']
no_edge_friction = int(options['no_edge_friction'])
lane_border_multiplier = int(options['lane_border_multiplier'])
min_tile_size = None
if options['min_tile_size']:
min_tile_size = float(options['min_tile_size'])
existing_cutlines = None
if options['existing_cutlines']:
existing_cutlines = options['existing_cutlines'].split(',')
tiles = options['output']
memory = int(options['memory'])
tiled = False
if options['tile_width']:
tiled = True
gscript.message(_("Using tiles processing for edge detection"))
width = int(options['tile_width'])
height = int(options['tile_height'])
overlap = int(options['overlap'])
processes = int(options['processes'])
global temp_maps
temp_maps = []
r = 'raster'
v = 'vector'
if existing_cutlines:
existingcutlinesmap = 'temp_icutlines_existingcutlinesmap_%i' % os.getpid()
if len(existing_cutlines) > 1:
gscript.run_command('v.patch',
input_=existing_cutlines,
output=existingcutlinesmap,
quiet=True,
overwrite=True)
existing_cutlines=existingcutlinesmap
gscript.run_command('v.to.rast',
input_=existing_cutlines,
output=existingcutlinesmap,
use='val',
type_='line,boundary',
overwrite=True,
quiet=True)
temp_maps.append([existingcutlinesmap, r])
temp_edge_map = "temp_icutlines_edgemap_%d" % os.getpid()
temp_maps.append([temp_edge_map, r])
gscript.message(_("Creating edge map"))
if edge_detection_algorithm == 'zc':
kwargs = {'input' : inputraster,
'output' : temp_edge_map,
'width_' : int(options['zc_width']),
'threshold' : float(options['zc_threshold']),
'quiet' : True}
if tiled:
grd = GridModule('i.zc',
width=width,
height=height,
overlap=overlap,
processes=processes,
split=False,
**kwargs)
grd.run()
else:
gscript.run_command('i.zc',
**kwargs)
elif edge_detection_algorithm == 'canny':
if not gscript.find_program('i.edge', '--help'):
message = _("You need to install the addon i.edge to use ")
message += _("the Canny edge detector.\n")
message += _(" You can install the addon with 'g.extension i.edge'")
gscript.fatal(message)
kwargs = {'input' : inputraster,
'output' : temp_edge_map,
'low_threshold' : float(options['canny_low_threshold']),
'high_threshold' : float(options['canny_high_threshold']),
'sigma' : float(options['canny_sigma']),
'quiet' : True}
if tiled:
grd = GridModule('i.edge',
width=width,
height=height,
overlap=overlap,
processes=processes,
split=False,
**kwargs)
grd.run()
else:
gscript.run_command('i.edge',
**kwargs)
else:
gscript.fatal("Only zero-crossing and Canny available as edge detection algorithms.")
region = gscript.region()
gscript.message(_("Finding cutlines in both directions"))
nsrange = float(region.n - region.s - region.nsres)
ewrange = float(region.e - region.w - region.ewres)
if nsrange > ewrange:
hnumber_lines = number_lines
vnumber_lines = int(number_lines * (nsrange / ewrange))
else:
vnumber_lines = number_lines
hnumber_lines = int(number_lines * (ewrange / nsrange))
# Create the lines in horizonal direction
nsstep = float(region.n - region.s - region.nsres) / hnumber_lines
hpointsy = [((region.n - i * nsstep) - region.nsres / 2.0) for i in range(0, hnumber_lines+1)]
hlanepointsy = [y - nsstep / 2.0 for y in hpointsy]
hstartpoints = listzip([region.w + 0.2 * region.ewres] * len(hpointsy), hpointsy)
hstoppoints = listzip([region.e - 0.2 * region.ewres] * len(hpointsy), hpointsy)
hlanestartpoints = listzip([region.w + 0.2 * region.ewres] * len(hlanepointsy), hlanepointsy)
hlanestoppoints = listzip([region.e - 0.2 * region.ewres] * len(hlanepointsy), hlanepointsy)
hlanemap = 'temp_icutlines_hlanemap_%i' % os.getpid()
temp_maps.append([hlanemap, v])
temp_maps.append([hlanemap, r])
os.environ['GRASS_VERBOSE'] = '0'
new = VectorTopo(hlanemap)
new.open('w')
for line in listzip(hlanestartpoints,hlanestoppoints):
new.write(geom.Line(line), cat=1)
new.close()
del os.environ['GRASS_VERBOSE']
gscript.run_command('v.to.rast',
input_=hlanemap,
output=hlanemap,
use='val',
type_='line',
overwrite=True,
quiet=True)
hbasemap = 'temp_icutlines_hbasemap_%i' % os.getpid()
temp_maps.append([hbasemap, r])
# Building the cost maps using the following logic
# - Any pixel not on an edge, nor on an existing cutline gets a
# no_edge_friction cost, or no_edge_friction_cost x 10 if there are
# existing cutlines
# - Any pixel on an edge gets a cost of 1 if there are no existing cutlines,
# and a cost of no_edge_friction if there are
# - A lane line gets a very high cost (lane_border_multiplier x cost of no
# edge pixel - the latter depending on the existence of cutlines).
mapcalc_expression = "%s = " % hbasemap
mapcalc_expression += "if(isnull(%s), " % hlanemap
if existing_cutlines:
mapcalc_expression += "if(%s == 0 && isnull(%s), " % (temp_edge_map, existingcutlinesmap)
mapcalc_expression += "%i, " % (no_edge_friction * 10)
mapcalc_expression += "if(isnull(%s), %s, 1))," % (existingcutlinesmap, no_edge_friction)
mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction * 10)
else:
mapcalc_expression += "if(%s == 0, " % temp_edge_map
mapcalc_expression += "%i, " % no_edge_friction
mapcalc_expression += "1), "
mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction)
gscript.run_command('r.mapcalc',
expression=mapcalc_expression,
quiet=True,
overwrite=True)
hcumcost = 'temp_icutlines_hcumcost_%i' % os.getpid()
temp_maps.append([hcumcost, r])
hdir = 'temp_icutlines_hdir_%i' % os.getpid()
temp_maps.append([hdir, r])
# Create the lines in vertical direction
ewstep = float(region.e - region.w - region.ewres) / vnumber_lines
vpointsx = [((region.e - i * ewstep) - region.ewres / 2.0) for i in range(0, vnumber_lines+1)]
vlanepointsx = [x + ewstep / 2.0 for x in vpointsx]
vstartpoints = listzip(vpointsx, [region.n - 0.2 * region.nsres] * len(vpointsx))
vstoppoints = listzip(vpointsx, [region.s + 0.2 * region.nsres] * len(vpointsx))
vlanestartpoints = listzip(vlanepointsx, [region.n - 0.2 * region.nsres] * len(vlanepointsx))
vlanestoppoints = listzip(vlanepointsx, [region.s + 0.2 * region.nsres] * len(vlanepointsx))
vlanemap = 'temp_icutlines_vlanemap_%i' % os.getpid()
temp_maps.append([vlanemap, v])
temp_maps.append([vlanemap, r])
os.environ['GRASS_VERBOSE'] = '0'
new = VectorTopo(vlanemap)
new.open('w')
for line in listzip(vlanestartpoints,vlanestoppoints):
new.write(geom.Line(line), cat=1)
new.close()
del os.environ['GRASS_VERBOSE']
gscript.run_command('v.to.rast',
input_=vlanemap,
output=vlanemap,
use='val',
type_='line',
overwrite=True,
quiet=True)
vbasemap = 'temp_icutlines_vbasemap_%i' % os.getpid()
temp_maps.append([vbasemap, r])
mapcalc_expression = "%s = " % vbasemap
mapcalc_expression += "if(isnull(%s), " % vlanemap
if existing_cutlines:
mapcalc_expression += "if(%s == 0 && isnull(%s), " % (temp_edge_map, existingcutlinesmap)
mapcalc_expression += "%i, " % (no_edge_friction * 10)
mapcalc_expression += "if(isnull(%s), %s, 1))," % (existingcutlinesmap, no_edge_friction)
mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction * 10)
else:
mapcalc_expression += "if(%s == 0, " % temp_edge_map
mapcalc_expression += "%i, " % no_edge_friction
mapcalc_expression += "1), "
mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction)
gscript.run_command('r.mapcalc',
expression=mapcalc_expression,
quiet=True,
overwrite=True)
vcumcost = 'temp_icutlines_vcumcost_%i' % os.getpid()
temp_maps.append([vcumcost, r])
vdir = 'temp_icutlines_vdir_%i' % os.getpid()
temp_maps.append([vdir, r])
if processes > 1:
pmemory = memory / 2.0
rcv = gscript.start_command('r.cost',
input_=vbasemap,
startcoordinates=vstartpoints,
stopcoordinates=vstoppoints,
output=vcumcost,
outdir=vdir,
memory=pmemory,
quiet=True,
overwrite=True)
rch = gscript.start_command('r.cost',
input_=hbasemap,
startcoordinates=hstartpoints,
stopcoordinates=hstoppoints,
output=hcumcost,
outdir=hdir,
memory=pmemory,
quiet=True,
overwrite=True)
rcv.wait()
rch.wait()
else:
gscript.run_command('r.cost',
input_=vbasemap,
startcoordinates=vstartpoints,
stopcoordinates=vstoppoints,
output=vcumcost,
outdir=vdir,
memory=memory,
quiet=True,
overwrite=True)
gscript.run_command('r.cost',
input_=hbasemap,
startcoordinates=hstartpoints,
stopcoordinates=hstoppoints,
output=hcumcost,
outdir=hdir,
memory=memory,
quiet=True,
overwrite=True)
hlines = 'temp_icutlines_hlines_%i' % os.getpid()
temp_maps.append([hlines, r])
vlines = 'temp_icutlines_vlines_%i' % os.getpid()
temp_maps.append([vlines, r])
if processes > 1:
rdh = gscript.start_command('r.drain',
input_=hcumcost,
direction=hdir,
startcoordinates=hstoppoints,
output=hlines,
flags='d',
quiet=True,
overwrite=True)
rdv = gscript.start_command('r.drain',
input_=vcumcost,
direction=vdir,
startcoordinates=vstoppoints,
output=vlines,
flags='d',
quiet=True,
overwrite=True)
rdh.wait()
rdv.wait()
else:
gscript.run_command('r.drain',
input_=hcumcost,
direction=hdir,
startcoordinates=hstoppoints,
output=hlines,
flags='d',
quiet=True,
overwrite=True)
gscript.run_command('r.drain',
input_=vcumcost,
direction=vdir,
startcoordinates=vstoppoints,
output=vlines,
flags='d',
quiet=True,
overwrite=True)
# Combine horizonal and vertical lines
temp_raster_tile_borders = 'temp_icutlines_raster_tile_borders_%i' % os.getpid()
temp_maps.append([temp_raster_tile_borders, r])
gscript.run_command('r.patch',
input_=[hlines,vlines],
output=temp_raster_tile_borders,
quiet=True,
overwrite=True)
gscript.message(_("Creating vector polygons"))
# Create vector polygons
# First we need to shrink the region a bit to make sure that all vector
# points / lines fall within the raster
gscript.use_temp_region()
gscript.run_command('g.region',
s=region.s+region.nsres,
e=region.e-region.ewres,
quiet=True)
region_map = 'temp_icutlines_region_map_%i' % os.getpid()
temp_maps.append([region_map, v])
temp_maps.append([region_map, r])
gscript.run_command('v.in.region',
output=region_map,
type_='line',
quiet=True,
overwrite=True)
gscript.del_temp_region()
gscript.run_command('v.to.rast',
input_=region_map,
output=region_map,
use='val',
type_='line',
quiet=True,
overwrite=True)
temp_raster_polygons = 'temp_icutlines_raster_polygons_%i' % os.getpid()
temp_maps.append([temp_raster_polygons, r])
gscript.run_command('r.patch',
input_=[temp_raster_tile_borders,region_map],
output=temp_raster_polygons,
quiet=True,
overwrite=True)
temp_raster_polygons_thin = 'temp_icutlines_raster_polygons_thin_%i' % os.getpid()
temp_maps.append([temp_raster_polygons_thin, r])
gscript.run_command('r.thin',
input_=temp_raster_polygons,
output=temp_raster_polygons_thin,
quiet=True,
overwrite=True)
# Create a series of temporary map names as we have to go
# through several steps until we reach the final map.
temp_vector_polygons1 = 'temp_icutlines_vector_polygons1_%i' % os.getpid()
temp_maps.append([temp_vector_polygons1, v])
temp_vector_polygons2 = 'temp_icutlines_vector_polygons2_%i' % os.getpid()
temp_maps.append([temp_vector_polygons2, v])
temp_vector_polygons3 = 'temp_icutlines_vector_polygons3_%i' % os.getpid()
temp_maps.append([temp_vector_polygons3, v])
temp_vector_polygons4 = 'temp_icutlines_vector_polygons4_%i' % os.getpid()
temp_maps.append([temp_vector_polygons4, v])
gscript.run_command('r.to.vect',
input_=temp_raster_polygons_thin,
output=temp_vector_polygons1,
type_='line',
flags='t',
quiet=True,
overwrite=True)
# Erase all category values from the lines
gscript.run_command('v.category',
input_=temp_vector_polygons1,
op='del',
cat='-1',
output=temp_vector_polygons2,
quiet=True,
overwrite=True)
# Transform lines to boundaries
gscript.run_command('v.type',
input_=temp_vector_polygons2,
from_type='line',
to_type='boundary',
output=temp_vector_polygons3,
quiet=True,
overwrite=True)
# Add centroids
gscript.run_command('v.centroids',
input_=temp_vector_polygons3,
output=temp_vector_polygons4,
quiet=True,
overwrite=True)
# If a threshold is given erase polygons that are too small
if min_tile_size:
gscript.run_command('v.clean',
input_=temp_vector_polygons4,
tool='rmarea',
threshold=min_tile_size,
output=tiles,
quiet=True,
overwrite=True)
else:
gscript.run_command('g.copy',
vect=[temp_vector_polygons4,tiles],
quiet=True,
overwrite=True)
gscript.vector_history(tiles)
if __name__ == "__main__":
options, flags = gscript.parser()
atexit.register(cleanup)
main()

38
opendm/grass/compute_cutline.grass
Wyświetl plik

@ -1,38 +0,0 @@
# orthophoto_file: input GeoTIFF raster file
# crop_area_file: input vector polygon file delimiting the safe area for processing
# number_lines: number of cutlines on the smallest side of the orthophoto for computing the final cutline
# max_concurrency: maximum number of parallel processes to use
# memory: maximum MB of memory to use
# ------
# output: If successful, prints the full path to the cutlines file. Otherwise it prints "error"
# Import orthophoto (green band only)
r.external band=2 input="${orthophoto_file}" output=ortho --overwrite
# Import crop area
v.in.ogr input="${crop_area_file}" output=crop_area --overwrite
g.region vector=crop_area
# Generate cutlines
i.cutlinesmod.py --overwrite input=ortho output=cutline number_lines=${number_lines} edge_detection=zc no_edge_friction=20 lane_border_multiplier=1000000 tile_width=1024 tile_height=1024 overlap=20 processes=${max_concurrency} memory=${memory}
#v.out.ogr input=cutline output="cutline_raw.gpkg" format=GPKG
# Select cutlines that are within crop area
v.select ainput=cutline binput=crop_area output=result operator=within
# Export
v.out.ogr input=result output="result.gpkg" format=GPKG --overwrite
# Merge all geometries, select only the largest one (remove islands)
ogr2ogr -f GPKG -overwrite -explodecollections -dialect SQLite -sql "SELECT ST_Union(geom) FROM result ORDER BY ST_AREA(geom) DESC LIMIT 1" cutline.gpkg result.gpkg
# Add new line output in case the last command didn't.
echo ""
if [ -e "cutline.gpkg" ]; then
echo "$$(pwd)/cutline.gpkg"
else
echo "error"
fi

143
opendm/grass_engine.py
Wyświetl plik

@ -1,143 +0,0 @@
import shutil
import tempfile
import subprocess
import os
import sys
import time
from opendm import log
from opendm import system
import locale
from string import Template
class GrassEngine:
def __init__(self):
self.grass_binary = system.which('grass7') or \
system.which('grass72') or \
system.which('grass74') or \
system.which('grass76') or \
shutil.which('grass78') or \
shutil.which('grass80')
if self.grass_binary is None:
log.ODM_WARNING("Could not find a GRASS 7 executable. GRASS scripts will not work.")
else:
log.ODM_INFO("Initializing GRASS engine using {}".format(self.grass_binary))
def create_context(self, serialized_context = {}):
if self.grass_binary is None: raise GrassEngineException("GRASS engine is unavailable")
return GrassContext(self.grass_binary, **serialized_context)
class GrassContext:
def __init__(self, grass_binary, tmpdir = None, template_args = {}, location = None, auto_cleanup=True):
self.grass_binary = grass_binary
if tmpdir is None:
tmpdir = tempfile.mkdtemp('_grass_engine')
self.tmpdir = tmpdir
self.template_args = template_args
self.location = location
self.auto_cleanup = auto_cleanup
def get_cwd(self):
return self.tmpdir
def add_file(self, filename, source, use_as_location=False):
param = os.path.splitext(filename)[0] # filename without extension
dst_path = os.path.abspath(os.path.join(self.get_cwd(), filename))
with open(dst_path, 'w') as f:
f.write(source)
self.template_args[param] = dst_path
if use_as_location:
self.set_location(self.template_args[param])
return dst_path
def add_param(self, param, value):
self.template_args[param] = value
def set_location(self, location):
"""
:param location: either a "epsg:XXXXX" string or a path to a geospatial file defining the location
"""
if not location.lower().startswith('epsg:'):
location = os.path.abspath(location)
self.location = location
def execute(self, script):
"""
:param script: path to .grass script
:return: script output
"""
if self.location is None: raise GrassEngineException("Location is not set")
script = os.path.abspath(script)
# Create grass script via template substitution
try:
with open(script) as f:
script_content = f.read()
except FileNotFoundError:
raise GrassEngineException("Script does not exist: {}".format(script))
tmpl = Template(script_content)
# Write script to disk
if not os.path.exists(self.get_cwd()):
os.mkdir(self.get_cwd())
with open(os.path.join(self.get_cwd(), 'script.sh'), 'w') as f:
f.write(tmpl.substitute(self.template_args))
# Execute it
log.ODM_INFO("Executing grass script from {}: {} --tmp-location {} --exec bash script.sh".format(self.get_cwd(), self.grass_binary, self.location))
env = os.environ.copy()
env["GRASS_ADDON_PATH"] = env.get("GRASS_ADDON_PATH", "") + os.path.abspath(os.path.join("opendm/grass/addons"))
env["LC_ALL"] = "C.UTF-8"
filename = os.path.join(self.get_cwd(), 'output.log')
with open(filename, 'wb') as writer, open(filename, 'rb', 1) as reader:
p = subprocess.Popen([self.grass_binary, '--tmp-location', self.location, '--exec', 'bash', 'script.sh'],
cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=writer, env=env)
while p.poll() is None:
sys.stdout.write(reader.read().decode('utf8'))
time.sleep(0.5)
# Read the remaining
sys.stdout.write(reader.read().decode('utf8'))
out, err = p.communicate()
out = out.decode('utf-8').strip()
if p.returncode == 0:
return out
else:
raise GrassEngineException("Could not execute GRASS script {} from {}: {}".format(script, self.get_cwd(), err))
def serialize(self):
return {
'tmpdir': self.tmpdir,
'template_args': self.template_args,
'location': self.location,
'auto_cleanup': self.auto_cleanup
}
def cleanup(self):
if os.path.exists(self.get_cwd()):
shutil.rmtree(self.get_cwd())
def __del__(self):
if self.auto_cleanup:
self.cleanup()
class GrassEngineException(Exception):
pass
def cleanup_grass_context(serialized_context):
ctx = grass.create_context(serialized_context)
ctx.cleanup()
grass = GrassEngine()

1
requirements.txt
Wyświetl plik

@ -28,3 +28,4 @@ scikit-image==0.17.2
scipy==1.5.4
xmltodict==0.12.0
fpdf2==2.2.0rc2
Shapely==1.7.1

2
snap/snapcraft.yaml
Wyświetl plik

@ -39,7 +39,6 @@ parts:
- gdal-bin
- gfortran # to build scipy
- git
- grass-core
- libboost-log-dev
- libgdal-dev
- libgeotiff-dev
@ -57,7 +56,6 @@ parts:
- swig3.0
stage-packages:
- gdal-bin
- grass-core
- libboost-log1.71.0
- libgdal26
- libgeotiff5

1
stages/odm_orthophoto.py
Wyświetl plik

@ -127,7 +127,6 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
bounds_file_path,
cutline_file,
args.max_concurrency,
tmpdir=os.path.join(tree.odm_orthophoto, "grass_cutline_tmpdir"),
scale=0.25)
orthophoto.compute_mask_raster(tree.odm_orthophoto_tif, cutline_file,