change TPI method

Author: yzh211

pyosp/_tpi.py

@@ -15,42 +15,30 @@ def __init__(self, point_coord, raster, radius):
         """
         self.p = point_coord
         self.raster = raster
-        self.geoTransform = self.raster.GetGeoTransform()
-        self.radius = int(radius / self.geoTransform[1])
-        
-        # change nodata to nan, assuming a small value was given 
-        self.rasterMatrix = self.raster.ReadAsArray()
-        if not np.min(self.rasterMatrix) > -1e10:
-            self.rasterMatrix[self.rasterMatrix < -1e10] = np.nan
+        self.cols = raster.RasterXSize
+        self.rows = raster.RasterYSize
+        self.geoTransform = raster.GetGeoTransform()
+        self.radiusInPixel = int(radius / self.geoTransform[1])
 
     def point_position(self):
         x = int((self.p[0] - self.geoTransform[0]) / self.geoTransform[1])
         y = int((self.geoTransform[3] - self.p[1]) / -self.geoTransform[5])
         return y, x
 
-    def raster_window(self):
+    def avg_window(self):
         py, px = self.point_position()
-        
-        #pad nan to range out of raster
-        rasterPad = np.pad(self.rasterMatrix, (self.radius+1,), 
-                           mode='constant', constant_values=(np.nan,))  
-                
-        return rasterPad[py:(py+2*self.radius+1), 
-                         px:(px+2*self.radius+1)]
+        xmin = max(0, px-self.radiusInPixel)
+        xmax = min(self.cols, px+self.radiusInPixel+1)
+        ymin = max(0, py-self.radiusInPixel)
+        ymax = min(self.rows, py+self.radiusInPixel+1)
+        arr = self.raster.ReadAsArray(xoff=xmin, yoff=ymin, xsize=xmax-xmin, ysize=ymax-ymin)
+        avg = (np.nansum(arr)-self.point_value()) / (arr.size-1)
+        return avg 
 
     def point_value(self):
-        return self.rasterMatrix[self.point_position()]
-    
-    @property
-    def index(self):
-        outer = self.window * self.raster_window()
-        inner = self.point_value()
-        return inner - (np.nansum(outer)/np.count_nonzero(~np.isnan(outer)))
+        py, px =self.point_position()
+        return self.raster.ReadAsArray(xoff=px, yoff=py, xsize=1, ysize=1)[0]
 
     @property
-    def window(self):
-        win = np.ones((2*self.radius+1, 2*self.radius+1))
-        r_y, r_x = win.shape[0]//2, win.shape[1]//2
-        win[r_y, r_x] = 0    # remove the central cell
-        return win
-    
+    def value(self):
+        return self.point_value() - self.avg_window() 

pyosp/cirsp/tpi_cir.py

@@ -63,7 +63,7 @@ def _radial_lines(self):
                 p = [self.center.x+dx, self.center.y+dy]
 
                 p_elev = Point_elevation(p, self.raster).value
-                p_tpi= Tpi(p, self.raster, self.tpi_radius).index
+                p_tpi= Tpi(p, self.raster, self.tpi_radius).value
 
                 if not (
                 (self.rasterXmin <= p[0] <= self.rasterXmax) and 

pyosp/curvsp/tpi_curv.py

@@ -5,6 +5,7 @@
 from ..util import pairwise, progressBar
 from .._tpi import Tpi
 from .base_curv import Base_curv
+import gdal
 
 class Tpi_curv(Base_curv):
     """TPI-based swath profile characterization.
@@ -29,13 +30,13 @@ class Tpi_curv(Base_curv):
     def __init__(self, line, raster, width, 
                  tpi_radius, min_tpi=float("-inf"), max_tpi=float("inf"),
                  line_stepsize=None, cross_stepsize=None):
-        self.tpi_radius = tpi_radius 
+        self.tpi_radius = tpi_radius
         self.min_tpi = min_tpi
         self.max_tpi = max_tpi
-        
+
         super(Tpi_curv,self).__init__(line, raster, width,
                                       line_stepsize, cross_stepsize)
-        
+
     def __repr__(self):
         return("{}".format(self.__class__.__name__))
 
@@ -104,8 +105,7 @@ def _swath_left(self, p1, p2, endPoint=False):
                 if hw >= self.width/2:
                     break
 
-            p_index = Tpi(p_left, self.raster, self.tpi_radius).index
-            
+            p_index = Tpi(p_left, self.raster, self.tpi_radius).value
             if self.min_tpi <= p_index <= self.max_tpi:
                 transect_temp.insert(0,p_left)
             else:
@@ -148,8 +148,7 @@ def _swath_right(self, p1, p2, endPoint=False):
                 if hw >= self.width/2:
                     break
 
-            p_index = Tpi(p_right, self.raster, self.tpi_radius).index                
-            
+            p_index = Tpi(p_right, self.raster, self.tpi_radius).value
             if self.min_tpi <= p_index <= self.max_tpi:
                 transect_temp.append(p_right)
             else:

pyosp/tests/__init__.py

@@ -0,0 +1,109 @@
+# -*- coding: utf-8 -*-
+
+import pytest
+import os, sys
+from collections import namedtuple
+from pyosp import *
+
+sys.path.append(os.path.dirname(os.path.abspath(__file__)))
+dat = os.path.join(os.path.dirname(os.path.abspath(__file__)), "../datasets/")
+
+homo_line = os.path.join(dat, 'homo_baseline.shp')
+homo_raster = os.path.join(dat, 'homo_mount.tif')
+
+cir_center = os.path.join(dat, 'center.shp')
+cir_raster = os.path.join(dat, 'crater.tif')
+
+@pytest.fixture(scope='module')
+def base_homo(**kwargs):
+    def _base_homo(**kwargs):
+        return Base_curv(line=homo_line, raster=homo_raster, width=100, **kwargs)
+    
+    return _base_homo
+
+@pytest.fixture(scope='module')
+def orig_homo(**kwargs):
+    def _orig_homo(**kwargs):
+        return Orig_curv(line=homo_line, raster=homo_raster, width=100, **kwargs)
+    
+    return _orig_homo
+
+@pytest.fixture(scope='module')
+def elev_homo(**kwargs):
+    def _elev_homo(**kwargs):
+        return Elev_curv(line=homo_line, raster=homo_raster, width=100,
+                         min_elev=0.01,
+                         **kwargs)
+    
+    return _elev_homo
+
+@pytest.fixture(scope='module')
+def slope_homo(**kwargs):
+    def _slope_homo(**kwargs):
+        return Slope_curv(line=homo_line, raster=homo_raster, width=100,
+                          min_slope=1,
+                          **kwargs)
+    
+    return _slope_homo
+
+@pytest.fixture(scope='module')
+def tpi_homo(**kwargs):
+    def _tpi_homo(**kwargs):
+        return Tpi_curv(line=homo_line, raster=homo_raster, width=100,
+                        tpi_radius=50, min_tpi=-5,
+                        **kwargs)
+    
+    return _tpi_homo
+
+@pytest.fixture(scope='module')
+def base_cir(**kwargs):
+    def _base_cir(**kwargs):
+        return Base_cir(cir_center, cir_raster, radius=80,
+                        ng_start=0, ng_end=300,
+                        ng_stepsize=1, radial_stepsize=None, **kwargs)
+    
+    return _base_cir
+
+@pytest.fixture(scope='module')
+def orig_cir(**kwargs):
+    def _orig_cir(**kwargs):
+        return Orig_cir(cir_center, cir_raster, radius=80,
+                        ng_start=0, ng_end=300,
+                        ng_stepsize=1, radial_stepsize=None, **kwargs)
+    
+    return _orig_cir
+
+@pytest.fixture(scope='module')
+def elev_cir(**kwargs):
+    def _elev_cir(**kwargs):
+        return Elev_cir(cir_center, cir_raster, radius=80,
+                        min_elev=4,
+                        ng_start=0, ng_end=300,
+                        ng_stepsize=1, radial_stepsize=None, **kwargs)
+    
+    return _elev_cir
+
+@pytest.fixture(scope='module')
+def slope_cir(**kwargs):
+    def _slope_cir(**kwargs):
+        return Slope_cir(cir_center, cir_raster, radius=80,
+                         min_slope=13,
+                         ng_start=0, ng_end=300,
+                         ng_stepsize=1, radial_stepsize=None, **kwargs)
+    
+    return _slope_cir
+
+@pytest.fixture(scope='module')
+def tpi_cir(**kwargs):
+    def _tpi_cir(**kwargs):
+        return Tpi_cir(cir_center, cir_raster, radius=80,
+                       tpi_radius=50, min_tpi=2,
+                       ng_start=0, ng_end=300,
+                       ng_stepsize=1, radial_stepsize=None, **kwargs)
+    
+    return _tpi_cir
+
+    
+
+    
+    

pyosp/tests/conftest.py

@@ -0,0 +1,15 @@
+# -*- coding: utf-8 -*-
+
+import os, sys
+
+sys.path.append(os.path.dirname(os.path.abspath(__file__)))
+dat = os.path.join(os.path.dirname(os.path.abspath(__file__)), "../datasets/")
+
+class TestBaseCir:
+
+    def test_initial(self, base_cir):
+        """Test the setup"""
+        base = base_cir()
+        assert base.radial_stepsize == base.raster.GetGeoTransform()[1]
+        assert len(base.distance) == base.radius // base.radial_stepsize + 1
+        

pyosp/tests/test_cirBase.py

@@ -0,0 +1,82 @@
+# -*- coding: utf-8 -*-
+
+import pytest
+import os, sys
+import pyosp
+import numpy as np
+
+sys.path.append(os.path.dirname(os.path.abspath(__file__)))
+dat = os.path.join(os.path.dirname(os.path.abspath(__file__)), "../datasets/")
+
+class TestCir:
+    def test_elev_cir(self, elev_cir):
+        """
+        Data INSIDE of border should fall in geo-parameter range.
+        OUTSIDE should not.
+        """
+        elev = elev_cir()
+        lines = elev.lines
+        
+        p_in = [x[-1] for x in lines[90:181]]
+        gradient = [np.radians(ng) for ng in range(90,181)]
+        p_dat_in = []
+        p_dat_out = []
+        for ind, point in enumerate(p_in):
+            p_out = [point[0] + np.cos(gradient[ind])*elev.radial_stepsize,
+                     point[1] + np.sin(gradient[ind])*elev.radial_stepsize]
+            
+            dat_in = pyosp.Point_elevation(point, elev.raster).value
+            dat_out = pyosp.Point_elevation(p_out, elev.raster).value
+            
+            p_dat_in.append(dat_in)
+            p_dat_out.append(dat_out)
+            
+        assert all(i < 4 for i in p_dat_out)
+        
+    def test_slope_cir(self, slope_cir):
+        """
+        Data INSIDE of border should fall in geo-parameter range.
+        OUTSIDE should not.
+        """
+        slope = slope_cir()
+        lines = slope.lines
+        cell_res = slope.raster.GetGeoTransform()[1]
+        
+        p_in = [x[-1] for x in lines[90:181]]
+        gradient = [np.radians(ng) for ng in range(90,181)]
+        p_dat_in = []
+        p_dat_out = []
+        for ind, point in enumerate(p_in):
+            p_out = [point[0] + np.cos(gradient[ind])*slope.radial_stepsize,
+                     point[1] + np.sin(gradient[ind])*slope.radial_stepsize]
+            
+            dat_in = pyosp.Geo_slope(point, slope.raster, cell_res).value
+            dat_out = pyosp.Geo_slope(p_out, slope.raster, cell_res).value
+            
+            p_dat_in.append(dat_in)
+            p_dat_out.append(dat_out)
+            
+        assert all(i < 13 for i in p_dat_out)
+        
+    def test_tpi_cir(self, tpi_cir):
+        """
+        Data INSIDE of border should fall in geo-parameter range.
+        OUTSIDE should not.
+        """
+        tpi = tpi_cir()
+        lines = tpi.lines
+        
+        p_in = [x[-1] for x in lines[90:181]]
+        gradient = [np.radians(ng) for ng in range(90,181)]
+        p_dat_in = []
+        p_dat_out = []
+        for ind, point in enumerate(p_in):
+            p_out = [point[0] + np.cos(gradient[ind])*tpi.radial_stepsize,
+                     point[1] + np.sin(gradient[ind])*tpi.radial_stepsize]
+            
+            dat_in = pyosp.Tpi(point, tpi.raster, tpi.tpi_radius).value
+            dat_out = pyosp.Tpi(p_out, tpi.raster, tpi.tpi_radius).value
+            p_dat_in.append(dat_in)
+            p_dat_out.append(dat_out)
+            
+        assert all(i < 2 for i in p_dat_out)

pyosp/tests/test_cirSP.py

@@ -0,0 +1,42 @@
+# -*- coding: utf-8 -*-
+
+import pytest
+import os, sys
+from pyosp import point_coords
+
+sys.path.append(os.path.dirname(os.path.abspath(__file__)))
+dat = os.path.join(os.path.dirname(os.path.abspath(__file__)), "../datasets/")
+
+class TestBaseHomo:
+    # @pytest.mark.parametrize('base_homo', [
+    #     dict(line_stepsize=None, cross_stepsize=None),
+    #     dict(line_stepsize=100)
+    #     ], indirect=True)
+    def test_initial(self, base_homo):
+        """Test the setup"""
+        base = base_homo(line_stepsize=None, cross_stepsize=None)
+        assert base.line_stepsize == base.cell_res
+        assert base.cross_stepsize == base.cell_res
+        assert len(base.distance) == base.line.length // base.line_stepsize + 1
+        assert len(base.line_p) == len(base.distance)
+        
+    def test_segment(self, base_homo):
+        """Test the start and end of chosen segment:
+            start == line.length / 2
+            end  == line.length / 4
+        """
+        start_end = os.path.join(dat, 'homo_start_end.shp')
+        coords = point_coords(start_end)
+        base = base_homo(line_stepsize=None, cross_stepsize=None)
+        start_ind, end_ind = base._segment(coords[1], coords[0])
+        assert abs(start_ind-len(base.distance)/4) <= 3
+        assert abs(end_ind-len(base.distance)/2) <= 3
+        
+        # if given values are distances
+        start_distance = base.distance[len(base.distance)//4]
+        end_distance = base.distance[len(base.distance)//2]
+        start_ind, end_ind = base._segment(start_distance, end_distance)
+        assert abs(start_ind-len(base.distance)/4) <= 3
+        assert abs(end_ind-len(base.distance)/2) <= 3
+        
+