Ensured that tile format and variable names consistent

biota/IO.py

@@ -89,8 +89,6 @@ def outputGeoTiff(data, filename, geo_t, proj, output_dir = os.getcwd(), dtype =
         ds.GetRasterBand(1).WriteArray(data)
     ds = None
 
-
-
 
 
 def _buildMap(fig, ax, data, lat, lon, title ='', cbartitle = '', vmin = 10., vmax = 60., cmap = 'YlGn'):
@@ -102,10 +100,10 @@ def _buildMap(fig, ax, data, lat, lon, title ='', cbartitle = '', vmin = 10., vm
 
     im = ax.imshow(data, vmin = vmin, vmax = vmax, cmap = cmap, interpolation = 'nearest')
 
-    ax.set_xticks(np.arange(0,4501,450))
-    ax.set_yticks(np.arange(0,4501,450))
-    ax.set_xticklabels(np.arange(lon, lon + 1.01, 0.1))
+    ax.set_yticks(np.arange(0,data.shape[0],data.shape[0]/10))
+    ax.set_xticks(np.arange(0,data.shape[1],data.shape[1]/10))
     ax.set_yticklabels(np.arange(lat, lat - 1.01, - 0.1))
+    ax.set_xticklabels(np.arange(lon, lon + 1.01, 0.1))
     ax.tick_params(labelsize = 5)
     ax.set_xlabel('Longitude', fontsize = 5)
     ax.set_ylabel('Latitude', fontsize = 5)
@@ -116,65 +114,60 @@ def _buildMap(fig, ax, data, lat, lon, title ='', cbartitle = '', vmin = 10., vm
     cbar.set_label(cbartitle, fontsize = 7)
 
 
-def overviewFigure(data_t1, data_t2, output_dir = os.getcwd(), output_name = 'overview'):
-    """overviewFigure(data_t1, data_t2, t1, t2, geo_t, output_dir = os.getcwd())
-    
-    Generate an overview image showing biomass and proportional biomass change for the tile being processed.
+def overviewFigure(tile_change, output = False, show = True):
+    """
+    Generate an overview image showing biomass and proportional biomass change for an ALOS tile.
     
     Args:
-        data_t1:
-        data_t2: 
-        output_name: Optionally specify an output string to precede output file. Defaults to 'overview'.
+        tile_change: An ALOS change object from biota.LoadChange()
     """
 
     import matplotlib.pyplot as plt
 
-    assert data_t1.geo_t == data_t2.geo_t, "The two ALOS tiles must be from the same location."
-
-    # Get upper left longitude and latitude from GeoMatrix
-    lon, lat = data_t1.lat, data_t1.lon
-
-    # Update masks to exclude areas outisde forest definition. Good for visualisation
-    AGB_t1 = data_t1.getAGB()
-    AGB_t2 = data_t2.getAGB()
+    # Load AGB, and update masks to exclude areas outisde forest definition. Good for visualisation
+    AGB_t1 = tile_change.data_t1.getAGB()
+    AGB_t2 = tile_change.data_t2.getAGB()
 
+    # Mask out areas < 10 tC/ha
     AGB_t1 = np.ma.array(AGB_t1, mask = np.logical_or(AGB_t1.mask, AGB_t1 < 10.))
     AGB_t2 = np.ma.array(AGB_t2, mask = np.logical_or(AGB_t2.mask, AGB_t1 < 10.))
 
-    AGB_change = (AGB_t2 - AGB_t1) / (data_t2.year - data_t1.year) # tC/ha/yr
+    AGB_change = AGB_t2 - AGB_t1
 
-    AGB_pcChange = 100 * (AGB_change / AGB_t1) # %/yr
+    AGB_pcChange = 100 * (AGB_change / AGB_t1) # %
 
     fig = plt.figure(figsize = (7, 6))
 
     # Plot a map of AGB at t1
     ax1 = fig.add_subplot(2, 2, 1)
-    _buildMap(fig, ax1, AGB_t1, lat, lon, title = 'AGB %s'%str(t1), cbartitle = 'tC/ha')
+    _buildMap(fig, ax1, AGB_t1, tile_change.lat, tile_change.lon, title = 'AGB %s'%str(tile_change.year_t1), cbartitle = 'tC/ha')
 
     # Plot a map of AGB at t2
     ax2 = fig.add_subplot(2, 2, 2)
-    _buildMap(fig, ax2, AGB_t2, lat, lon, title = 'AGB %s'%str(t2), cbartitle = 'tC/ha')    
+    _buildMap(fig, ax2, AGB_t2, tile_change.lat, tile_change.lon, title = 'AGB %s'%str(tile_change.year_t2), cbartitle = 'tC/ha')    
 
     # Plot a map of absolute AGB change   
     ax3 = fig.add_subplot(2, 2, 3)
-    _buildMap(fig, ax3, AGB_change, lat, lon, title = 'AGB change (%s-%s)'%(str(t1),str(t2)),
-              cbartitle = 'tC/ha/yr', vmin = -10., vmax = 10., cmap = 'RdBu')    
+    _buildMap(fig, ax3, AGB_change, tile_change.lat, tile_change.lon, title = 'AGB change (%s-%s)'%(str(tile_change.data_t1.year),str(tile_change.year_t2)),
+              cbartitle = 'tC/ha', vmin = -10., vmax = 10., cmap = 'RdBu')    
 
     # Plot a map of % AGB change
     ax4 = fig.add_subplot(2, 2, 4)
-    _buildMap(fig, ax4, AGB_pcChange, lat, lon, title = 'AGB change (%s-%s)'%(str(t1),str(t2)),
-              cbartitle = '%/yr', vmin = -50., vmax = 50., cmap = 'RdBu')    
+    _buildMap(fig, ax4, AGB_pcChange, tile_change.lat, tile_change.lon, title = 'AGB change (%s-%s)'%(str(tile_change.data_t1.year),str(tile_change.year_t2)),
+              cbartitle = '%', vmin = -50., vmax = 50., cmap = 'RdBu')    
 
     plt.tight_layout()
 
-    # Determine filename
-    hem_NS = 'S' if lat < 0 else 'N'
-    hem_EW = 'W' if lon < 0 else 'E'
+    if output:
+        output_pattern = tile_change.output_pattern.replace('.tif','.png')
 
-    output_path = '%s/%s_%s%s.png'%(output_dir, output_name, data_t1.hem_NS + str(abs(lat)).zfill(2), 
-                                    data_t1.hem_EW + str(abs(lon)).zfill(3))
+        output_path = os.path.abspath(os.path.expanduser('%s/%s'%(tile_change.output_dir, output_pattern%('OverviewFigure'))))
 
-    plt.savefig(output_path, dpi = 150)
+        plt.savefig(output_path, dpi = 150)
+
+    if show:
+        plt.show()
+
     plt.close()
 
 

biota/calibrate.py

@@ -11,25 +11,27 @@
 
 import pdb
 
-def extractGamma0(shp, plot_field, agb_field, data_example, verbose = False, units = 'natural'):
+
+
+def extractGamma0(shp, plot_field, agb_field, tile_example, verbose = False, units = 'natural'):
     '''
     Extract gamma0 from ALOS tiles.
     
     Args:
         shp: A shapefile containing plot data
         plot_field: The shapefile field containing plot names
         agb_field: The shapefile field containing AGB estimates
-        data_example: An exemplar ALOS tile (from LoadTile()) which contains information on the processing chain.
+        tile_example: An exemplar ALOS tile (from LoadTile()) which contains information on the processing chain.
         #TODO: year?
     Returns:
         A dictionary containing plot names, AGB and gamma0 values
     '''
 
     # Use example data to get processing steps
-    downsample_factor = data_example.downsample_factor
-    lee_filter = data_example.lee_filter
-    year = data_example.year
-    dataloc = data_example.dataloc
+    downsample_factor = tile_example.downsample_factor
+    lee_filter = tile_example.lee_filter
+    year = tile_example.year
+    dataloc = tile_example.dataloc
 
     # Extract relevant info from shapefile
     plot_names = biota.mask.getField(shp, plot_field) 
@@ -54,21 +56,21 @@ def extractGamma0(shp, plot_field, agb_field, data_example, verbose = False, uni
         if verbose: print 'Doing lat: %s, lon: %s'%(str(lat), str(lon))
 
         # Load tile
-        data = biota.LoadTile(dataloc, lat, lon, year, downsample_factor = downsample_factor, lee_filter = lee_filter)
+        tile = biota.LoadTile(dataloc, lat, lon, year, downsample_factor = downsample_factor, lee_filter = lee_filter)
 
         # Get backscatter (both polarisations) and DOY
-        data_gamma0_hv = data.getGamma0(polarisation = 'HV', units = units)
-        data_gamma0_hh = data.getGamma0(polarisation = 'HH', units = units)
-        data_doy = data.getDOY()
+        data_gamma0_hv = tile.getGamma0(polarisation = 'HV', units = units)
+        data_gamma0_hh = tile.getGamma0(polarisation = 'HH', units = units)
+        data_doy = tile.getDOY()
 
         # Mask out each plot
-        plot_mask = biota.mask.rasterizeShapefile(data, shp, location_id = True, buffer_size = 0.)
+        plot_mask = biota.mask.rasterizeShapefile(tile, shp, location_id = True, buffer_size = 0.)
 
         # Extract values for each plot
         for n in np.unique(plot_mask[plot_mask != 0]):
 
-            # Get mask for plot and data
-            this_mask = np.logical_and(plot_mask == n, data.mask == False)
+            # Get mask for plot and tile
+            this_mask = np.logical_and(plot_mask == n, tile.mask == False)
 
             # Add metrics to output array
             gamma0_hv[n-1] = np.mean(data_gamma0_hv[this_mask])

biota/core.py

@@ -361,7 +361,7 @@ def updateMask(self, filename, buffer_size = 0., classes = []):
         if file_type == 'shp':
 
             # Rasterize the shapefile, optionally with a buffer
-            mask = biota.mask.rasterizeShapefile(self, filename, buffer_size = buffer_size)
+            mask = biota.mask.maskShapefile(self, filename, buffer_size = buffer_size)
 
         if file_type in ['tif', 'tiff']:
 

biota/mask.py

@@ -164,12 +164,12 @@ def getBBox(shp, field, value):
     return bbox
 
 
-def maskRaster(data, tif, classes = [], buffer_size = 0.):
+def maskRaster(tile, tif, classes = [], buffer_size = 0.):
     '''
     Extract a mask from a GeoTiff based on specified classes
     
     Args:
-        data: An ALOS object
+        tile: An ALOS tile (biota.LoadTile())
         tif: A GeoTiff file, with integer values
         classes: A list of values to add to be masked
         buffer_size: Optionally specify a buffer to add around maked pixels, in meters.
@@ -188,12 +188,12 @@ def maskRaster(data, tif, classes = [], buffer_size = 0.):
 
     # Create output file matching ALOS tile
     gdal_driver = gdal.GetDriverByName('MEM')
-    ds_dest = gdal_driver.Create('', data.ySize, data.xSize, 1, 3)
-    ds_dest.SetGeoTransform(data.geo_t)
-    ds_dest.SetProjection(data.proj)
+    ds_dest = gdal_driver.Create('', tile.ySize, tile.xSize, 1, 3)
+    ds_dest.SetGeoTransform(tile.geo_t)
+    ds_dest.SetProjection(tile.proj)
 
     # Reproject input GeoTiff to match the ALOS tile
-    gdal.ReprojectImage(ds_source, ds_dest, proj_source, data.proj, gdal.GRA_NearestNeighbour)
+    gdal.ReprojectImage(ds_source, ds_dest, proj_source, tile.proj, gdal.GRA_NearestNeighbour)
 
     # Load resampled image into memory
     tif_resampled = ds_dest.GetRasterBand(1).ReadAsArray()
@@ -204,20 +204,20 @@ def maskRaster(data, tif, classes = [], buffer_size = 0.):
     if buffer_size > 0.:
 
         # Determine the size of the buffer in pixels
-        buffer_px = int(buffer_size / ((data.xRes + data.yRes) / 2.))
+        buffer_px = int(buffer_size / ((tile.xRes + tile.yRes) / 2.))
 
         # Dilate the mask
         mask = dilateMask(mask, buffer_px)
 
     return mask
 
 
-def rasterizeShapefile(data, shp, buffer_size = 0., field = None, value = None, location_id = False):
+def maskShapefile(tile, shp, buffer_size = 0., field = None, value = None, location_id = False):
     """
     Rasterize points, lines or polygons from a shapefile to match ALOS mosaic data.
         
     Args:
-        data: An ALOS object
+        tile: An ALOS tile (biota.LoadTile())
         shp: Path to a shapefile consisting of points, lines and/or polygons. This does not have to be in the same projection as ds
         buffer_size: Optionally specify a buffer to add around features of the shapefile, in meters.
         field: Optionally specify a single shapefile field to extract (you must also specify its value)
@@ -232,19 +232,15 @@ def rasterizeShapefile(data, shp, buffer_size = 0., field = None, value = None,
     from osgeo import gdalnumeric
 
     assert np.logical_or(np.logical_and(field == None, value == None), np.logical_and(field != None, value != None)), "If specifying field or value, both must be defined. At present, field = %s and value = %s"%(str(field), str(value))
-
-    # Determine size of buffer to place around lines/polygons
-    #buffer_px = int(round(buffer_size / data.geo_t[1]))
-
+
     # Determine the size of the buffer in degrees
-    buffer_size_degrees = buffer_size / (((data.xRes * data.xSize) + (data.yRes * data.ySize)) / 2.)
+    buffer_size_degrees = buffer_size / (((tile.xRes * tile.xSize) + (tile.yRes * tile.ySize)) / 2.)
 
     # Determine size of buffer to place around lines/polygons
-    #buffer_px = int(round(buffer_size / ((data.xRes + data.yRes) / 2.)))
-    buffer_px = int(round(buffer_size_degrees / data.geo_t[1]))
+    buffer_px = int(round(buffer_size_degrees / tile.geo_t[1]))
 
     # Create output image. Add a buffer around the image array equal to the maxiumum dilation size. This means that features just outside ALOS tile extent can contribute to dilated mask.
-    rasterPoly = Image.new("I", (data.ySize + (buffer_px * 2), data.xSize + (buffer_px * 2)), 0)
+    rasterPoly = Image.new("I", (tile.ySize + (buffer_px * 2), tile.xSize + (buffer_px * 2)), 0)
     rasterize = ImageDraw.Draw(rasterPoly)
 
     # The shapefile may not have the same CRS as ALOS mosaic data, so this will generate a function to reproject points.
@@ -272,10 +268,10 @@ def rasterizeShapefile(data, shp, buffer_size = 0., field = None, value = None,
         sxmax, symax, z = coordTransform.TransformPoint(sxmax, symax)
 
         # Go to the next record if out of bounds
-        if sxmax < data.geo_t[0] - buffer_size_degrees: continue
-        if sxmin > data.geo_t[0] + (data.geo_t[1] * data.ySize) + buffer_size_degrees: continue
-        if symax < data.geo_t[3] + (data.geo_t[5] * data.xSize) + buffer_size_degrees: continue
-        if symin > data.geo_t[3] - buffer_size_degrees: continue
+        if sxmax < tile.geo_t[0] - buffer_size_degrees: continue
+        if sxmin > tile.geo_t[0] + (tile.geo_t[1] * tile.ySize) + buffer_size_degrees: continue
+        if symax < tile.geo_t[3] + (tile.geo_t[5] * tile.xSize) + buffer_size_degrees: continue
+        if symin > tile.geo_t[3] - buffer_size_degrees: continue
 
         #Separate polygons with list indices
         n_parts = len(shape.parts) #Number of parts
@@ -297,7 +293,7 @@ def rasterizeShapefile(data, shp, buffer_size = 0., field = None, value = None,
                 lon, lat, z = coordTransform.TransformPoint(p[0], p[1])
 
                 # Then convert map to pixel coordinates using geo transform
-                pixels.append(_world2Pixel(data.geo_t, lon, lat, buffer_size = buffer_size_degrees))
+                pixels.append(_world2Pixel(tile.geo_t, lon, lat, buffer_size = buffer_size_degrees))
 
             # Draw the mask for this shape...
             # if a point...