Commit 7a7821f2 authored by Gaetano Raffaele's avatar Gaetano Raffaele
Browse files

ENH/FIX: major changes in roi-based workflow.

parent 1b46a606
......@@ -277,7 +277,7 @@ def getRasterExtent(ras_ds):
return ext
def getRasterExtentAsGeometry(ras_ds):
def getRasterExtentAsGeometry(ras_ds, to_wgs84=False):
""" getRasterExtentAsGeometry(ras_ds)
Gets the extent of a raster GDAL dataset as polygon (ogr.Geometry).
......@@ -294,8 +294,20 @@ def getRasterExtentAsGeometry(ras_ds):
ring.AddPoint(extent[0], extent[2])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
if to_wgs84:
ras_srs = osr.SpatialReference(wkt=ras_ds.GetProjection())
out_srs = osr.SpatialReference()
out_srs.ImportFromEPSG(4326)
poly.AssignSpatialReference(ras_srs)
poly.TransformTo(out_srs)
return poly
def getRasterEPSG(fn):
ds = gdal.Open(fn)
srs = osr.SpatialReference(wkt=ds.GetProjection())
return srs.GetAttrValue('AUTHORITY',1)
def getRasterExtentAsShapefile(ras, shp):
""" getRasterExtentAsShapefile(ras,shp)
......
#import argparse
import glob, os, sys, warnings
#from path import Path
#from os.path import join, isdir
import subprocess
from itertools import groupby
import gdal, ogr, osr
#import geopandas as gpd
import otbApplication as otb
from mtdUtils import queuedProcess, getBufferedGeographicExtent, getRasterExtentAsGeometry, getGeoJSONExtent_WGS84, randomword
from mtdUtils import getBufferedGeographicExtent, getRasterEPSG, getGeoJSONExtent_WGS84, getRasterExtentAsGeometry
import xml.etree.ElementTree as ET
#import shutil
import multiprocessing as mp
from path import Path
......@@ -19,9 +15,6 @@ import zipfile
from shutil import unpack_archive
import elevation
#from shapely.geometry import box
#import rasterio as rio
# TODO:
# - rasterize vector on
# - check possibilities to update filter outcore to make filtering on a temporal moving window
......@@ -119,7 +112,7 @@ def download_geoid(cache_dir=Path('~/.moringa/cache').expanduser(), verbose=True
return geoid_file
def S1CalibOrtho(fld, epsg, demfld, geoid, overwrite=False):
def S1CalibOrtho(fld, epsg=None, demfld=None, geoid=None, overwrite=False):
lst = glob.glob(fld + os.sep + 'measurement' + os.sep + '*.tiff')
sarcal = []
ortho = []
......@@ -145,6 +138,53 @@ def S1CalibOrtho(fld, epsg, demfld, geoid, overwrite=False):
return out
def S1CalibOrthoWithROI(infile, outfile, roi, demfld, geoid, lut='sigma', overwrite=False, ram=1024):
if checkRoiInS1Acquisition(infile, roi):
pipe = []
pipe.append(otb.Registry.CreateApplication('SARCalibration'))
pipe[-1].SetParameterString('in', infile)
pipe[-1].SetParameterString('lut', lut)
pipe[-1].SetParameterString('ram', str(ram))
pipe[-1].Execute()
pipe.append(otb.Registry.CreateApplication('OrthoRectification'))
pipe[-1].SetParameterInputImage('io.in', pipe[-2].GetParameterOutputImage('out'))
pipe[-1].SetParameterString('elev.dem', demfld)
pipe[-1].SetParameterString('elev.geoid', geoid)
pipe[-1].SetParameterString('map', 'epsg')
pipe[-1].SetParameterString('map.epsg.code', getRasterEPSG(roi))
pipe[-1].SetParameterInt('opt.gridspacing', 40)
pipe[-1].SetParameterString('opt.ram', str(ram))
pipe[-1].Execute()
pipe.append(otb.Registry.CreateApplication('Superimpose'))
pipe[-1].SetParameterInputImage('inm', pipe[-2].GetParameterOutputImage('io.out'))
pipe[-1].SetParameterString('inr', roi)
pipe[-1].SetParameterString('ram', str(ram))
pipe[-1].SetParameterString('out', outfile)
if (os.path.exists(outfile) and not overwrite):
print('Calibration and Orthorectification process skipped, file already exists:\n\t{}'.format(outfile))
else:
pipe[-1].ExecuteAndWriteOutput()
return outfile
else:
print('ROI not overlapping with S1 acquisition :\n\t{}'.format(infile))
return None
def checkRoiInS1Acquisition(s1file, roifile):
isIntersecting = False
dsv = gdal.Open(roifile)
roi_geom = getRasterExtentAsGeometry(dsv, to_wgs84=True)
dsi = ogr.Open(os.path.dirname(s1file) + os.sep + os.pardir + os.sep + 'preview' + os.sep + 'map-overlay.kml')
lyi = dsi.GetLayer()
for g in lyi:
in_geom = g.GetGeometryRef()
isIntersecting = in_geom.Intersects(roi_geom)
dsv = None
dsi = None
return isIntersecting
def preClip(infile, outfile, roi_buffered_vector, overwrite=False, verbose=True, ram=1024):
if Path(outfile).exists() and (not overwrite):
......@@ -185,15 +225,16 @@ def preClip(infile, outfile, roi_buffered_vector, overwrite=False, verbose=True,
# subprocess.check_call(cmd, shell=False)
def S1Calibration(infile, outfile, overwrite=False, verbose=True, ram=1024):
def S1Calibration(infile, outfile, lut='sigma', overwrite=False, verbose=True, ram=1024):
if Path(outfile).exists() and (not overwrite):
if verbose:
print('Calibration process skipped, file already exists:\n\t{}'.format(outfile))
return outfile
print(lut)
sarcal = otb.Registry.CreateApplication('SARCalibration')
sarcal.SetParameterString('in', infile)
sarcal.SetParameterString('lut', 'sigma')
sarcal.SetParameterString('lut', lut)
sarcal.SetParameterString('out', outfile)
sarcal.SetParameterString('ram', str(ram))
sarcal.ExecuteAndWriteOutput()
......@@ -604,7 +645,7 @@ def raster_epsg(x):
return int(proj.GetAttrValue('AUTHORITY', 1))
def s1process(indir, outdir, epsg=None,
roi=None, buffer=100, dem=None, geoid=None, orthofit=False,
roi=None, buffer=100, dem=None, geoid=None, orthofit=False, lut='sigma',
direction='descending', satellite='both', step='all', force_clip=False, features=None, skip_db_conversion=False,
features_by_date=False, calib_dir=None, cache_dir=Path('~/.moringa/cache').expanduser(),
ram=1024, ncores=1, overwrite=False, verbose=True):
......@@ -685,10 +726,8 @@ def s1process(indir, outdir, epsg=None,
if orthofit:
all_steps = ['calib', 'ortho', 'mosaic', 'filter']
if force_clip:
all_steps = ['clip'] + all_steps
else:
all_steps = ['clip', 'calib', 'ortho', 'mosaic', 'superimpose', 'filter']
all_steps = ['clip', 'calib', 'ortho', 'mosaic', 'filter']
if isinstance(step, str):
step = [step]
......@@ -770,8 +809,39 @@ def s1process(indir, outdir, epsg=None,
#raster_bbox(roi, buf=buffer, outfile=roi_buffered_vector)
getBufferedGeographicExtent(roi, buf=buffer, toFile=roi_buffered_vector, drvname='GeoJSON')
#### clipping
if (roi is not None) and ('clip' in step):
if (roi is not None and 'clip' in step and 'calib' in step and 'ortho' in step):
#### pipelined calibration + ortho + superimpose
ortho_dir.mkdir_p()
if dem_dir is None:
dem_dir = outdir / 'dem'
dem_dir.mkdir_p()
srtm_file = dem_dir / 'srtm.tiff'
download_srtm1(roi_buffered_vector, srtm_file, cache_dir)
if geoid is None:
geoid = download_geoid(cache_dir, verbose)
outfiles = [ortho_dir / infile.name.replace('.tiff', '_clipped_cal' + lut + '_ortho.tiff') for infile in proc_files]
mp_args = [(infile, outfile, roi, dem_dir, geoid, lut, overwrite, ram) for infile, outfile in
zip(proc_files, outfiles)]
"""
files = []
for arg in mp_args:
try:
f = S1CalibOrthoWithROI(*arg)
if f is not None:
files.append(f)
except Exception as e:
print(e)
print('Clip went wrong, file skipped: {}'.format(f))
"""
with mp.Pool(ncores) as pool:
files = pool.starmap(S1CalibOrthoWithROI, mp_args) # 1h for 8 tiles in parallel (8 threads)
proc_files = [x for x in files if x is not None]
else:
if (roi is not None and 'clip' in step):
### Clipping
if verbose:
print('Clipping process...')
clip_dir.mkdir_p()
......@@ -796,14 +866,15 @@ def s1process(indir, outdir, epsg=None,
print(e)
print('Clip went wrong, file skipped: {}'.format(f))
proc_files = files
#### calibration
if 'calib' in step:
if verbose:
print('Calibration process...')
calib_dir.mkdir_p()
outfiles = [calib_dir / infile.name.replace('.tiff', '_calsigma.tiff') for infile in proc_files]
mp_args = [(infile, outfile, overwrite, verbose, ram) for infile, outfile in zip(proc_files, outfiles)]
outfiles = [calib_dir / infile.name.replace('.tiff', '_cal' + lut + '.tiff') for infile in proc_files]
mp_args = [(infile, outfile, lut, overwrite, verbose, ram) for infile, outfile in zip(proc_files, outfiles)]
files = []
for arg in mp_args:
f = S1Calibration(*arg) # 6s/S2tile using 32threads
......@@ -831,7 +902,6 @@ def s1process(indir, outdir, epsg=None,
geoid = download_geoid(cache_dir, verbose)
outfiles = [ortho_dir / infile.name.replace('.tiff', '_ortho.tiff') for infile in proc_files]
print(geoid)
mp_args = [(infile, outfile, epsg, dem_dir, geoid, roi, orthofit, overwrite, verbose, ram) for infile, outfile in
zip(proc_files, outfiles)]
# files = []
......@@ -887,9 +957,10 @@ def parse_cmd_line(argv=sys.argv[1:]):
parser.add_argument("--orthofit", action='store_true', help="If activated it will use mode orthofit with ROI in OTB OrthoRectification app."
"Otherwise image is clipped before calibration step and "
"aligned (OTB Superimpose) on ROI grid after mosaic")
parser.add_argument("--force_clip", action='store_true', help="When orthofit mode is used, default is to calibrate the whole S1 image before"
"orthorectification. If set, the S1 scene will be pre-clipped before orhto"
"(introduces some distortion to the calibration process).")
parser.add_argument("--lut",
help='Specify lookup table to use for calibration (sigma, gamma, beta).',
choices=['sigma', 'gamma', 'beta'],
default='sigma')
parser.add_argument("--direction",
help="Direction of acquisition, either ascending, descending",
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment