RivaMap experimental scripts

pairwise_analysis.py

# -*- coding: utf-8 -*-
"""
Created on Sun Jun 12 22:11:50 2016

@author: leo

Computes the difference in channelness/islandness between two images
"""

import cv2
import numpy as np
from rivamap import singularity_index, georef

nrScales = 14
minScale = 1.2

im1path = '/home/leo/rivamap/data/waxlake/waxlake2015/ae68719ed2bc1560c66741cda9db6d61.nd.tif'
im2path = '/home/leo/rivamap/data/waxlake/waxlake2005/550158ebaa8e58141fcabb838f6b84c6.nd.tif'

I1 = cv2.imread(im1path, cv2.IMREAD_UNCHANGED)
gm1 = georef.loadGeoMetadata(im1path)

I2 = cv2.imread(im2path, cv2.IMREAD_UNCHANGED)
gm2 = georef.loadGeoMetadata(im2path)

filters = singularity_index.SingularityIndexFilters(minScale, nrScales)
psi1, widthMap1, orient1 = singularity_index.applyMMSI(I1, filters, togglePolarity = True)
psi2, widthMap2, orient2 = singularity_index.applyMMSI(I2, filters, togglePolarity = True)

psi_diff = psi1 - psi2

georef.saveAsGeoTiff(gm1, psi_diff, "island_diff_geotagged.TIF")

temporal_analysis.py

"""
Computes island and channel instability over time given a set of geotiff images
"""

import cv2
import numpy as np
from rivamap import singularity_index, georef
import glob, os

# All tiff files should be in their folders under the base directory
base_dir = "/home/leo/rivamap/data/ganges/"
mndwi_images = glob.glob(base_dir + '*/*.tif')

#year_range = range(1984, 2016)

#initialize psi_array
I1 = cv2.imread(mndwi_images[0], cv2.IMREAD_UNCHANGED)
psi_array = np.zeros((len(mndwi_images), I1.shape[0], I1.shape[1]), dtype='float32')
gm = georef.loadGeoMetadata(mndwi_images[0])

# island response change
filters = singularity_index.SingularityIndexFilters(nrScales=16)
for i in range(0, len(mndwi_images)):
    I1 = cv2.imread(mndwi_images[i], cv2.IMREAD_UNCHANGED)
    psi, widthMap, orient = singularity_index.applyMMSI(I1, filters, togglePolarity = True)
    psi_array[i,:,:] = psi

island_var = np.var(psi_array, axis=0)

# channel response change
filters = singularity_index.SingularityIndexFilters(nrScales=16)
for i in range(0, len(mndwi_images)):
    I1 = cv2.imread(mndwi_images[i], cv2.IMREAD_UNCHANGED)
    psi, widthMap, orient = singularity_index.applyMMSI(I1, filters, togglePolarity = False)
    psi_array[i,:,:] = psi

channel_var = np.var(psi_array, axis=0)

var_diff = channel_var - island_var

georef.saveAsGeoTiff(gm, var_diff, os.path.join(base_dir, "var_diff_geotagged.TIF"))