

def derive_datasets(self, input_dataset_dict, stack_output_info, tile_type_info):
    " Overrides abstract function in stacker class. Called in Stacker.stack_derived() function. \n        Creates PQA-masked NDVI stack\n        \n        Arguments:\n            nbar_dataset_dict: Dict keyed by processing level (e.g. ORTHO, NBAR, PQA, DEM)\n                containing all tile info which can be used within the function\n                A sample is shown below (including superfluous band-specific information):\n                \n{\n'NBAR': {'band_name': 'Visible Blue',\n    'band_tag': 'B10',\n    'end_datetime': datetime.datetime(2000, 2, 9, 23, 46, 36, 722217),\n    'end_row': 77,\n    'level_name': 'NBAR',\n    'nodata_value': -999L,\n    'path': 91,\n    'satellite_tag': 'LS7',\n    'sensor_name': 'ETM+',\n    'start_datetime': datetime.datetime(2000, 2, 9, 23, 46, 12, 722217),\n    'start_row': 77,\n    'tile_layer': 1,\n    'tile_pathname': '/g/data/v10/datacube/EPSG4326_1deg_0.00025pixel/LS7_ETM/150_-025/2000/LS7_ETM_NBAR_150_-025_2000-02-09T23-46-12.722217.tif',\n    'x_index': 150,\n    'y_index': -25},\n'ORTHO': {'band_name': 'Thermal Infrared (Low Gain)',\n     'band_tag': 'B61',\n     'end_datetime': datetime.datetime(2000, 2, 9, 23, 46, 36, 722217),\n     'end_row': 77,\n     'level_name': 'ORTHO',\n     'nodata_value': 0L,\n     'path': 91,\n     'satellite_tag': 'LS7',\n     'sensor_name': 'ETM+',\n     'start_datetime': datetime.datetime(2000, 2, 9, 23, 46, 12, 722217),\n     'start_row': 77,\n     'tile_layer': 1,\n     'tile_pathname': '/g/data/v10/datacube/EPSG4326_1deg_0.00025pixel/LS7_ETM/150_-025/2000/LS7_ETM_ORTHO_150_-025_2000-02-09T23-46-12.722217.tif',\n     'x_index': 150,\n     'y_index': -25},\n'PQA': {'band_name': 'Pixel Quality Assurance',\n    'band_tag': 'PQA',\n    'end_datetime': datetime.datetime(2000, 2, 9, 23, 46, 36, 722217),\n    'end_row': 77,\n    'level_name': 'PQA',\n    'nodata_value': None,\n    'path': 91,\n    'satellite_tag': 'LS7',\n    'sensor_name': 'ETM+',\n    'start_datetime': datetime.datetime(2000, 2, 9, 23, 46, 12, 722217),\n    'start_row': 77,\n    'tile_layer': 1,\n    'tile_pathname': '/g/data/v10/datacube/EPSG4326_1deg_0.00025pixel/LS7_ETM/150_-025/2000/LS7_ETM_PQA_150_-025_2000-02-09T23-46-12.722217.tif,\n    'x_index': 150,\n    'y_index': -25}\n}                \n                \n        Arguments (Cont'd):\n            stack_output_info: dict containing stack output information. \n                Obtained from stacker object. \n                A sample is shown below\n                \nstack_output_info = {'x_index': 144, \n                      'y_index': -36,\n                      'stack_output_dir': '/g/data/v10/tmp/ndvi',\n                      'start_datetime': None, # Datetime object or None\n                      'end_datetime': None, # Datetime object or None \n                      'satellite': None, # String or None \n                      'sensor': None} # String or None \n                      \n        Arguments (Cont'd):\n            tile_type_info: dict containing tile type information. \n                Obtained from stacker object (e.g: stacker.tile_type_dict[tile_type_id]). \n                A sample is shown below\n                \n{'crs': 'EPSG:4326',\n    'file_extension': '.tif',\n    'file_format': 'GTiff',\n    'format_options': 'COMPRESS=LZW,BIGTIFF=YES',\n    'tile_directory': 'EPSG4326_1deg_0.00025pixel',\n    'tile_type_id': 1L,\n    'tile_type_name': 'Unprojected WGS84 1-degree at 4000 pixels/degree',\n    'unit': 'degree',\n    'x_origin': 0.0,\n    'x_pixel_size': Decimal('0.00025000000000000000'),\n    'x_pixels': 4000L,\n    'x_size': 1.0,\n    'y_origin': 0.0,\n    'y_pixel_size': Decimal('0.00025000000000000000'),\n    'y_pixels': 4000L,\n    'y_size': 1.0}\n                            \n        Function must create one or more GDAL-supported output datasets. Useful functions in the\n        Stacker class include Stacker.get_pqa_mask(), but it is left to the coder to produce exactly\n        what is required for a single slice of the temporal stack of derived quantities.\n            \n        Returns:\n            output_dataset_info: Dict keyed by stack filename\n                containing metadata info for GDAL-supported output datasets created by this function.\n                Note that the key(s) will be used as the output filename for the VRT temporal stack\n                and each dataset created must contain only a single band. An example is as follows:\n{'/g/data/v10/tmp/ndvi/NDVI_stack_150_-025.vrt': \n    {'band_name': 'Normalised Differential Vegetation Index with PQA applied',\n    'band_tag': 'NDVI',\n    'end_datetime': datetime.datetime(2000, 2, 9, 23, 46, 36, 722217),\n    'end_row': 77,\n    'level_name': 'NDVI',\n    'nodata_value': None,\n    'path': 91,\n    'satellite_tag': 'LS7',\n    'sensor_name': 'ETM+',\n    'start_datetime': datetime.datetime(2000, 2, 9, 23, 46, 12, 722217),\n    'start_row': 77,\n    'tile_layer': 1,\n    'tile_pathname': '/g/data/v10/tmp/ndvi/LS7_ETM_NDVI_150_-025_2000-02-09T23-46-12.722217.tif',\n    'x_index': 150,\n    'y_index': -25}\n}\n        "
    log_multiline(logger.debug, input_dataset_dict, 'input_dataset_dict', '\t')
    rgb_bands = ('SWIR1', 'NIR', 'G')
    rgb_minmax = ((780, 5100), (200, 4500), (100, 2300))
    nbar_dataset_info = input_dataset_dict.get('NBAR')
    if (not nbar_dataset_info):
        log_multiline(logger.warning, input_dataset_dict, 'NBAR dict does not exist', '\t')
        return None
    lookup = BandLookup(data_cube=self, lookup_scheme_name='LANDSAT-LS5/7', tile_type_id=tile_type_info['tile_type_id'], satellite_tag=nbar_dataset_info['satellite_tag'], sensor_name=nbar_dataset_info['sensor_name'], level_name=nbar_dataset_info['level_name'])
    nbar_dataset_path = nbar_dataset_info['tile_pathname']
    output_tile_path = os.path.join(self.output_dir, re.sub('\\.\\w+$', '_RGB.tif', os.path.basename(nbar_dataset_path)))
    if os.path.exists(output_tile_path):
        logger.info('Skipping existing file %s', output_tile_path)
        return None
    if (not self.lock_object(output_tile_path)):
        logger.info('Skipping locked file %s', output_tile_path)
        return None
    input_dataset = gdal.Open(nbar_dataset_path)
    assert input_dataset, ('Unable to open dataset %s' % nbar_dataset_path)
    if ((nbar_dataset_info['satellite_tag'] == 'LS8') and (nbar_dataset_info['sensor_name'] == 'OLI')):
        pqa_mask = numpy.ones(shape=(input_dataset.RasterYSize, input_dataset.RasterXSize), dtype=numpy.bool)
        logger.debug('Work-around for LS8-OLI PQA issue applied: EVERYTHING PASSED')
    else:
        if (input_dataset_dict.get('PQA') is None):
            return
        pqa_mask = self.get_pqa_mask(pqa_dataset_path=input_dataset_dict['PQA']['tile_pathname'])
    log_multiline(logger.debug, pqa_mask, 'pqa_mask', '\t')
    gdal_driver = gdal.GetDriverByName('GTiff')
    output_dataset = gdal_driver.Create(output_tile_path, input_dataset.RasterXSize, input_dataset.RasterYSize, 3, gdal.GDT_Byte, ['INTERLEAVE=PIXEL', 'COMPRESS=LZW'])
    assert output_dataset, ('Unable to open output dataset %s' % output_dataset)
    output_dataset.SetGeoTransform(input_dataset.GetGeoTransform())
    output_dataset.SetProjection(input_dataset.GetProjection())
    for band_index in range(3):
        logger.debug('Processing %s band in layer %s as band %s', rgb_bands[band_index], lookup.band_no[rgb_bands[band_index]], (band_index + 1))
        scale = ((rgb_minmax[band_index][1] - rgb_minmax[band_index][0]) / 254.0)
        offset = (1.0 - (rgb_minmax[band_index][0] / scale))
        input_array = input_dataset.GetRasterBand(lookup.band_no[rgb_bands[band_index]]).ReadAsArray()
        log_multiline(logger.debug, input_array, 'input_array', '\t')
        output_array = ((input_array / scale) + offset).astype(numpy.byte)
        output_array[(input_array < rgb_minmax[band_index][0])] = 1
        output_array[(input_array > rgb_minmax[band_index][1])] = 255
        output_array[(~ pqa_mask)] = 0
        log_multiline(logger.debug, output_array, 'output_array', '\t')
        output_band = output_dataset.GetRasterBand((band_index + 1))
        output_band.WriteArray(output_array)
        output_band.SetNoDataValue(0)
        output_band.FlushCache()
    output_dataset.FlushCache()
    self.unlock_object(output_tile_path)
    logger.info('Finished writing RGB file %s', output_tile_path)
    return None

dataset/ETHPy150Open GeoscienceAustralia/agdc/examples/rgb_stacker.py IndexStacker.derive_datasets/original