From 0b075e1dedbbbc4d52410f646c93baa8fbaeb031 Mon Sep 17 00:00:00 2001
From: "raffaele.gaetano" <raffaele.gaetano@cirad.fr>
Date: Tue, 20 Nov 2018 15:08:50 +0100
Subject: [PATCH] FIX: Workaround for L8 south hemisphere / north projection
issue.
---
coregister.py | 19 ++++++++++++-------
1 file changed, 12 insertions(+), 7 deletions(-)
diff --git a/coregister.py b/coregister.py
index 5a63770..d3751aa 100644
--- a/coregister.py
+++ b/coregister.py
@@ -30,11 +30,11 @@ def main(argv):
src = args[1]
ds = gdal.Open(args[1])
- #prj = ds.GetProjection()
+ prj = ds.GetProjection()
gt = ds.GetGeoTransform()
- #srs = osr.SpatialReference()
- #srs.ImportFromWkt(prj)
- #code = srs.GetAuthorityCode(None)
+ srs = osr.SpatialReference()
+ srs.ImportFromWkt(prj)
+ code = srs.GetAuthorityCode(None)
gsp = str(int(2 * round(max(abs(gt[1]), abs(gt[5])))))
ds = None
@@ -105,7 +105,7 @@ def main(argv):
homo_sample_rate = float(maxgcp)/float(npts)
homoGeo2Pixel(tmpfld + 'homo.txt', args[1], tmpfld + 'homocorr.txt', sample_rate=homo_sample_rate)
- cmd = ['otbcli_GenerateRPCSensorModel', '-inpoints', tmpfld + 'homocorr.txt', '-outgeom', tmphomo]
+ cmd = ['otbcli_GenerateRPCSensorModel', '-inpoints', tmpfld + 'homocorr.txt', '-outgeom', tmphomo, '-map', 'epsg', 'map.epsg.code', '4326']
subprocess.call(cmd, shell=sh)
# Remove Cartographic Information from source (temporary)
@@ -116,7 +116,12 @@ def main(argv):
ds.SetProjection('')
ds = None
- cmd = ['otbcli_OrthoRectification', '-io.in \"' + args[1] + '?&geom=' + tmphomo + '\"', '-io.out', tmpfld + 'sits_band_coreg.tif', 'uint16', '-map', 'epsg', '-map.epsg.code', '4326', '-opt.gridspacing', gsp]
+ tmp_code = code
+ # Workaround for southern hemisphere images with NORTH projection
+ if tmp_code.startswith('326') and gt[3] < 0:
+ warnings.warn("Original image is in southern hemisphere with north UTM projection. Using south projection for temporary output.")
+ tmp_code = tmp_code[0:2] + '7' + tmp_code[3:]
+ cmd = ['otbcli_OrthoRectification', '-io.in \"' + args[1] + '?&geom=' + tmphomo + '\"', '-io.out', tmpfld + 'sits_band_coreg.tif', 'uint16', '-map', 'epsg', '-map.epsg.code', tmp_code, '-opt.gridspacing', gsp]
cmdex = ' '.join(cmd)
subprocess.call(cmdex, shell=True)
@@ -144,7 +149,7 @@ def main(argv):
ds.SetProjection('')
ds = None
- cmd = ['otbcli_OrthoRectification', '-io.in \"' + f + '?&geom=' + tmphomo + '\"', '-io.out', tmpfld + fn[0] + '_coreg.tif', 'uint8', '-map', 'epsg', '-map.epsg.code', '4326', '-opt.gridspacing', gsp, '-interpolator nn']
+ cmd = ['otbcli_OrthoRectification', '-io.in \"' + f + '?&geom=' + tmphomo + '\"', '-io.out', tmpfld + fn[0] + '_coreg.tif', 'uint8', '-map', 'epsg', '-map.epsg.code', tmp_code, '-opt.gridspacing', gsp, '-interpolator nn']
cmdex = ' '.join(cmd)
subprocess.call(cmdex, shell=True)
--
GitLab