diff --git a/coregister.py b/coregister.py index 5a63770491f8eb9b7aa77ba3b97c0ecc1b682b07..d3751aa982c8d7394224362d21eee78e4bb3693f 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)