FIX: Workaround for L8 south hemisphere / north projection issue.

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)