ffs

Author: Ciaran1981

geospatial_learn/handyplots.py

@@ -1,16 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-Created on Wed May  6 14:21:38 2015
-@author: Ciaran Robb
-author: Ciaran Robb
-Research Associate in Earth Observation
-Centre for Landscape and Climate Research (CLCR)
-Department of Geography, University of Leicester, University Road, Leicester, 
-LE1 7RH, UK 
-
-If you use code to publish work cite/acknowledge me and authors of libs as 
-appropriate 
-
+The handy plots module - these plots may be/ may have been handy
 
 
 """
@@ -61,7 +51,7 @@ def plot_classif_report(trueVals, predVals, labels=None, target_names=None,
 
     cbVl = dF.values
     mn = np.round(cbVl.min(), decimals=2)
-    mx= np.round(cbVl.max(), decimals=2)
+    mx = np.round(cbVl.max(), decimals=2)
     del cbVl
 
 

geospatial_learn/learning.py

@@ -211,7 +211,7 @@ def objective(trial, X, y, cv, group, score=scr):
         param_grid = {
             #"device_type": trial.suggest_categorical("device_type", ['gpu']),
            #'metric': 'rmse', 
-            'random_state': 48,
+            'random_state': 42,
             'n_estimators': 20000,
             'reg_alpha': trial.suggest_loguniform('reg_alpha', 1e-3, 10.0),
             'reg_lambda': trial.suggest_loguniform('reg_lambda', 1e-3, 10.0),

geospatial_learn/raster.py

@@ -1523,7 +1523,7 @@ def combine_scene(scl, c_scn, blocksize = 256):
     inDataset.FlushCache()
     inDataset = None
 
-def polygonize(inRas, outPoly, outField=None,  mask = True, band = 1, 
+def polygonize(inRas, outPoly, outField=None,  mask=True, band=1, 
                filetype="ESRI Shapefile"):
 
     """ 
@@ -1585,8 +1585,8 @@ def polygonize(inRas, outPoly, outField=None,  mask = True, band = 1,
 
     if outField is None:
         dst_fieldname = 'DN'
-        fd = ogr.FieldDefn( dst_fieldname, ogr.OFTInteger )
-        dst_layer.CreateField( fd )
+        fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
+        dst_layer.CreateField(fd)
         dst_field = dst_layer.GetLayerDefn().GetFieldIndex(dst_fieldname)
 
 
@@ -1767,7 +1767,8 @@ def set_bandnames(inras, names):
     rds.FlushCache()
     rds = None
 
-def rasterize(inShp, inRas, outRas, field=None, fmt="Gtiff"):
+def rasterize(inShp, inRas, outRas, field=None, fmt="Gtiff", 
+              dtype=gdal.GDT_Int32):
 
     """ 
     Rasterize a polygon to the extent & geo transform of another raster
@@ -1787,6 +1788,10 @@ def rasterize(inShp, inRas, outRas, field=None, fmt="Gtiff"):
     
     fmt: the gdal image format
     
+    dtype: int
+            gdal dtype for output raster - consider the dtytpe of the field
+            you are burning
+    
     """
 
 
@@ -1796,7 +1801,7 @@ def rasterize(inShp, inRas, outRas, field=None, fmt="Gtiff"):
     # the usual 
 
     outDataset = _copy_dataset_config(inDataset, FMT=fmt, outMap=outRas,
-                         dtype = gdal.GDT_Int32, bands=1)
+                         dtype=dtype, bands=1)
 
 
     vds = ogr.Open(inShp)
@@ -1806,6 +1811,7 @@ def rasterize(inShp, inRas, outRas, field=None, fmt="Gtiff"):
     if field == None:
         gdal.RasterizeLayer(outDataset, [1], lyr, burn_values=[1])
     else:
+        print('rasterizing', field)
         gdal.RasterizeLayer(outDataset, [1], lyr, options=["ATTRIBUTE="+field])
 
     outDataset.FlushCache()
@@ -2022,8 +2028,8 @@ def clip_raster(inRas, inShp, outRas, cutline=False, fmt='GTiff'):
         rds1 = None
 
 
-def fill_nodata(inRas, maxSearchDist=5, smoothingIterations=1, 
-                bands=[1]):
+def fill_nodata(inRas, maxSearchDist=3, smoothingIterations=0, 
+                bands=[1], nodata=0):
 
     """
     fill no data using gdal
@@ -2043,12 +2049,16 @@ def fill_nodata(inRas, maxSearchDist=5, smoothingIterations=1,
     bands: list of ints
             the bands to process      
     
+    nodata: int/float
+            the nodata value (it may not be set)  
+    
     """
 
     rds = gdal.Open(inRas, gdal.GA_Update)
 
     for band in tqdm(bands):
         bnd = rds.GetRasterBand(band)
+        bnd.SetNoDataValue(nodata)
         gdal.FillNodata(targetBand=bnd, maskBand=None, 
                          maxSearchDist=maxSearchDist, 
                          smoothingIterations=smoothingIterations)
@@ -2525,7 +2535,8 @@ def _copy_dataset_config(inDataset, FMT = 'Gtiff', outMap = 'copy',
         x_pixels,
         y_pixels,
         bands,
-        dtype)
+        dtype,
+        options=['COMPRESS=LZW'])
 
     outDataset.SetGeoTransform((
         x_min,    # 0

geospatial_learn/shape.py

@@ -1025,8 +1025,7 @@ def zonal_stats(inShp, inRas, band, bandname, layer=None, stat = 'mean',
         elif stat == 'min':
             feature_stats = float(masked.min())
         elif stat == 'mean':
-            feature_stats = float(masked.mean())
-            
+            feature_stats = float(masked.mean())     
         elif stat == 'max':
             feature_stats = float(masked.max())
         elif stat == 'median':
@@ -1038,9 +1037,20 @@ def zonal_stats(inShp, inRas, band, bandname, layer=None, stat = 'mean',
         elif stat == 'count':
             feature_stats = int(masked.count())
         elif stat == 'perc':
-            total = masked.shape[0]* masked.shape[1]
-            perc = masked.count() / total
-            feature_stats = int(np.round(perc*100))  
+            # this would be the perc of the bbox
+            #total = masked.shape[0]* masked.shape[1]
+            #perc = masked.count() / total
+            # mask does not help here
+            #TODO crap hack needs sorted
+            srcn = np.count_nonzero(src_array)
+            rvcn =  np.count_nonzero(rv_array)
+            if srcn > rvcn:
+                perc = 1
+            elif srcn == 0:
+                perc = 0
+            else:
+                perc = srcn / rvcn
+            feature_stats = float(np.round(perc, decimals=2))  
         elif stat == 'var':
             feature_stats = float(masked.var())
         elif stat == 'skew':
@@ -1265,7 +1275,7 @@ def zonal_frac(inShp, inRas, band, bandname, layer=None,
             unique = src_array[0]
             count = np.array([1])
         else:
-            unique, count = np.unique(src_array[rv_array>0], return_counts=True)
+            unique, count = np.unique(src_array[rv_array>0], return_counts=True,)
 
         entry = [feat.GetFID(), unique.tolist(), count.tolist()]
 
@@ -2550,16 +2560,18 @@ def zonal_point(inShp, inRas, field, band=1, nodata_value=0, write_stat=True):
     vds = ogr.Open(inShp, 1)  # TODO maybe open update if we want to write stats
     vlyr = vds.GetLayer(0)
 
-    if write_stat != None:
-        # if the field exists leave it as ogr is a pain with dropping it
-        # plus can break the file
-        if _fieldexist(vlyr, field) == False:
-            vlyr.CreateField(ogr.FieldDefn(field, ogr.OFTReal))
+    # writing via OGR is glacial
+    # if write_stat != None:
+    #     # if the field exists leave it as ogr is a pain with dropping it
+    #     # plus can break the file
+    #     if _fieldexist(vlyr, field) == False:
+    #         vlyr.CreateField(ogr.FieldDefn(field, ogr.OFTReal))
 
 
 
     feat = vlyr.GetNextFeature()
     features = np.arange(vlyr.GetFeatureCount())
+    outdata = []
 
     for label in tqdm(features):
 
@@ -2584,19 +2596,21 @@ def zonal_point(inShp, inRas, field, band=1, nodata_value=0, write_stat=True):
                 # unlikely but if none will have no data in the attribute table
                 continue
             outval =  float(src_array.max())
-            
-#            if write_stat != None:
-            feat.SetField(field, outval)
-            vlyr.SetFeature(feat)
+            outdata.append(outval)
+# #            if write_stat != None:
+#             feat.SetField(field, outval)
+#             vlyr.SetFeature(feat)
             feat = vlyr.GetNextFeature()
-        
+    
+    gdf = gpd.read_file(inShp)
+    gdf[field] = np.array(outdata)
     if write_stat != None:
-        vlyr.SyncToDisk()
-
-
-
+        gdf.to_file(inShp)
+    #     vlyr.SyncToDisk()
     vds = None
     rds = None
+    
+    return gdf
 
 def zonal_point_stk(inShp, inRas, fields, write_stat=True):
 
@@ -2661,7 +2675,6 @@ def zonal_point_stk(inShp, inRas, fields, write_stat=True):
             mx, my = geom.GetX(), geom.GetY()  
 
             # Convert from map to pixel coordinates.
-            # No rotation but for this that should not matter
             px = int((mx - rgt[0]) / rgt[1])
             py = int((my - rgt[3]) / rgt[5])