Python raster.Raster类代码示例

 def test_qgis_raster_layer_loading(self):
     """Test that reading from QgsRasterLayer works."""
     # This line is the cause of the problem:
     qgis_layer = QgsRasterLayer(RASTER_BASE + '.tif', 'test')
     layer = Raster(data=qgis_layer)
     qgis_extent = qgis_layer.dataProvider().extent()
     qgis_extent = [qgis_extent.xMinimum(), qgis_extent.yMinimum(),
                    qgis_extent.xMaximum(), qgis_extent.yMaximum()]
     layer_exent = layer.get_bounding_box()
     self.assertListEqual(
         layer_exent, qgis_extent,
         'Expected %s extent, got %s' % (qgis_extent, layer_exent))

    def test_convert_to_qgis_raster_layer(self):
        """Test that converting to QgsVectorLayer works."""
        if qgis_imported:
            # Create vector layer
            keywords = read_keywords(RASTER_BASE + '.keywords')
            layer = Raster(data=RASTER_BASE + '.tif', keywords=keywords)

            # Convert to QgsRasterLayer
            qgis_layer = layer.as_qgis_native()
            qgis_extent = qgis_layer.dataProvider().extent()
            qgis_extent = [qgis_extent.xMinimum(), qgis_extent.yMinimum(),
                           qgis_extent.xMaximum(), qgis_extent.yMaximum()]
            layer_exent = layer.get_bounding_box()
            self.assertListEqual(
                layer_exent, qgis_extent,
                'Expected %s extent, got %s' % (qgis_extent, layer_exent))

def write_raster_data(data, projection, geotransform, filename, keywords=None):
    """Write array to raster file with specified metadata and one data layer

    Input:
        data: Numpy array containing grid data
        projection: WKT projection information
        geotransform: 6 digit vector
                      (top left x, w-e pixel resolution, rotation,
                       top left y, rotation, n-s pixel resolution).
                       See e.g. http://www.gdal.org/gdal_tutorial.html
        filename: Output filename
        keywords: Optional dictionary

    Note: The only format implemented is GTiff and the extension must be .tif
    """

    R = Raster(data, projection, geotransform, keywords=keywords)
    R.write_to_file(filename)

#.........这里部分代码省略.........
                    (grid_point - covered_exposure_top_left) / (
                        covered_exposure_dimension)).astype(int)
                new_covered_exposure_data[index[1]][index[0]] = 0

        # Estimate number of people in need of evacuation
        if self.use_affected_field:
            affected_population = tr(
                'People within hazard field ("%s") of value "%s"') % (
                    self.hazard_class_attribute,
                    ','.join([
                        unicode(hazard_class) for
                        hazard_class in self.hazard_class_mapping[self.wet]
                    ]))
        else:
            affected_population = tr('People within any hazard polygon.')

        self.affected_population[affected_population] = (
            total_affected_population)

        self.total_population = int(
            numpy.nansum(self.exposure.layer.get_data(scaling=False)))
        self.unaffected_population = (
            self.total_population - self.total_affected_population)

        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]

        # Create style
        colours = ['#FFFFFF', '#38A800', '#79C900', '#CEED00',
                   '#FFCC00', '#FF6600', '#FF0000', '#7A0000']
        classes = create_classes(
            new_covered_exposure_data.flat[:], len(colours))

        # check for zero impact
        if total_affected_population == 0:
            message = no_population_impact_message(self.question)
            raise ZeroImpactException(message)

        interval_classes = humanize_class(classes)
        # Define style info for output polygons showing population counts
        style_classes = []
        for i in xrange(len(colours)):
            style_class = dict()
            style_class['label'] = create_label(interval_classes[i])
            if i == 1:
                label = create_label(
                    interval_classes[i],
                    tr('Low Population [%i people/cell]' % classes[i]))
            elif i == 4:
                label = create_label(
                    interval_classes[i],
                    tr('Medium Population [%i people/cell]' % classes[i]))
            elif i == 7:
                label = create_label(
                    interval_classes[i],
                    tr('High Population [%i people/cell]' % classes[i]))
            else:
                label = create_label(interval_classes[i])

            style_class['label'] = label
            style_class['quantity'] = classes[i]
            style_class['colour'] = colours[i]
            style_class['transparency'] = 0
            style_classes.append(style_class)

        # Override style info with new classes and name
        style_info = dict(
            target_field=None,
            style_classes=style_classes,
            style_type='rasterStyle')

        impact_data = self.generate_data()

        extra_keywords = {
            'target_field': self.target_field,
            'map_title': self.metadata().key('map_title'),
            'legend_notes': self.metadata().key('legend_notes'),
            'legend_units': self.metadata().key('legend_units'),
            'legend_title': self.metadata().key('legend_title'),
            'affected_population': total_affected_population,
            'total_population': self.total_population,
            'total_needs': self.total_needs
        }

        impact_layer_keywords = self.generate_impact_keywords(extra_keywords)

        # Create raster layer and return
        impact_layer = Raster(
            data=new_covered_exposure_data,
            projection=covered_exposure.get_projection(),
            geotransform=covered_exposure.get_geotransform(),
            name=self.metadata().key('layer_name'),
            keywords=impact_layer_keywords,
            style_info=style_info)

        impact_layer.impact_data = impact_data
        self._impact = impact_layer
        return impact_layer

def start(west,north,east,south, since, until=None, data_dir=None, population=None):
    
    bbox = (west, north, east, south)

    year, month, day = [int(x) for x in since.split('-')]
    since = datetime.date(year, month, day)

    if not isinstance(until, datetime.date):
        year, month, day = [int(x) for x in until.split('-')]
        until = datetime.date(year, month, day)
    else:
        until = until

    # Make sure the inputs are divisible by 10.
    for item in bbox:
        msg = "%d is not divisible by 10." % item
        assert int(item) % 10 == 0, msg

    the_viewports = viewports(bbox)
    the_timespan = timespan(since, until)

    data_dir = os.path.abspath(data_dir)

    if not os.path.exists(data_dir):
        os.mkdir(data_dir)

    print 'Downloading layers per day'
    # Download the layers for the given viewport and timespan.
    download(the_viewports, the_timespan, data_dir)

    print 'Merging layers per day'
    merged_files = merge(the_timespan, data_dir)

    flood_filename = os.path.join(data_dir, 'flood_severity.tif')

    if not os.path.exists(flood_filename):
        if len(merged_files) > 0:
            # Add all the pixels with a value higher than 3.
            #accumulate(merged_files, flood_filename, threshold=3)
            flooded = _flood_severity(merged_files)
            flooded.write_to_file(flood_filename)

            subprocess.call(['gdal_merge.py',
                     '-co', 'compress=packbits',
                     '-o', 'flood_severity_compressed.tif',
                     '-ot', 'Byte',
                     flood_filename], stdout=open(os.devnull, 'w'))
            os.remove(flood_filename)
            os.rename('flood_severity_compressed.tif', flood_filename)
        else:
            raise Exception('No merged files found for %s' % the_timespan)
    
    population_file = os.path.join(data_dir, population)
    population_object = Raster(population_file)
    # get population bbox
    pop_bbox = population_object.get_bounding_box()

    # get resolutions and pick the best
    pop_resolution = population_object.get_resolution()[0]

    hazard_object = Raster(flood_filename)
    hazard_resolution = hazard_object.get_resolution()[0]
    hazard_bbox = hazard_object.get_bounding_box()

    if pop_bbox[0] > bbox[0] and pop_bbox[1] > bbox[1] and pop_bbox[2] < bbox[2] and pop_bbox[3] < bbox[3]:
        hazard_file = clip(flood_filename, pop_bbox, cellSize=pop_resolution)
        exposure_layer = population_file
    else:
        hazard_file = clip(flood_filename, hazard_bbox, cellSize=pop_resolution)
        exposure_layer = clip(population_file, hazard_bbox, cellSize=None)    

    basename, ext = os.path.splitext(hazard_file)
    keywords_file = basename + '.keywords'

    if not os.path.exists(keywords_file):
        with open(keywords_file, 'w') as f:
            f.write(FLOOD_KEYWORDS)

    impact = calculate(hazard_file, exposure_layer)

    impact.write_to_file('impact.tif')

    count = impact.keywords['count']
    pretty_date = until.strftime('%a %d, %b %Y')
    print pretty_date, "|", "People affected: %s / %s" % (count, impact.keywords['total'])

#.........这里部分代码省略.........
        exposure_data = self.exposure.layer.get_data(nan=True, scaling=True)
        if has_no_data(exposure_data):
            self.no_data_warning = True

        # Make 3 data for each zone. Get the value of the exposure if the
        # exposure is in the hazard zone, else just assign 0
        low_exposure = numpy.where(hazard_data < low_t, exposure_data, 0)
        medium_exposure = numpy.where(
            (hazard_data >= low_t) & (hazard_data < medium_t),
            exposure_data, 0)
        high_exposure = numpy.where(
            (hazard_data >= medium_t) & (hazard_data <= high_t),
            exposure_data, 0)
        impacted_exposure = low_exposure + medium_exposure + high_exposure

        # Count totals
        self.total_population = int(numpy.nansum(exposure_data))
        self.affected_population[
            tr('Population in high hazard areas')] = int(
                numpy.nansum(high_exposure))
        self.affected_population[
            tr('Population in medium hazard areas')] = int(
                numpy.nansum(medium_exposure))
        self.affected_population[
            tr('Population in low hazard areas')] = int(
                numpy.nansum(low_exposure))
        self.unaffected_population = (
            self.total_population - self.total_affected_population)

        # check for zero impact
        if self.total_affected_population == 0:
            message = no_population_impact_message(self.question)
            raise ZeroImpactException(message)

        # Don't show digits less than a 1000
        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]
        total_needs = self.total_needs

        # Style for impact layer
        colours = [
            '#FFFFFF', '#38A800', '#79C900', '#CEED00',
            '#FFCC00', '#FF6600', '#FF0000', '#7A0000']
        classes = create_classes(impacted_exposure.flat[:], len(colours))
        interval_classes = humanize_class(classes)
        style_classes = []

        for i in xrange(len(colours)):
            style_class = dict()
            if i == 1:
                label = create_label(
                    interval_classes[i],
                    tr('Low Population [%i people/cell]' % classes[i]))
            elif i == 4:
                label = create_label(
                    interval_classes[i],
                    tr('Medium Population [%i people/cell]' % classes[i]))
            elif i == 7:
                label = create_label(
                    interval_classes[i],
                    tr('High Population [%i people/cell]' % classes[i]))
            else:
                label = create_label(interval_classes[i])
            style_class['label'] = label
            style_class['quantity'] = classes[i]
            style_class['transparency'] = 0
            style_class['colour'] = colours[i]
            style_classes.append(style_class)

        style_info = dict(
            target_field=None,
            style_classes=style_classes,
            style_type='rasterStyle')

        impact_data = self.generate_data()

        extra_keywords = {
            'map_title': self.metadata().key('map_title'),
            'legend_notes': self.metadata().key('legend_notes'),
            'legend_units': self.metadata().key('legend_units'),
            'legend_title': self.metadata().key('legend_title'),
            'total_needs': total_needs
        }

        impact_layer_keywords = self.generate_impact_keywords(extra_keywords)

        # Create raster object and return
        impact_layer = Raster(
            data=impacted_exposure,
            projection=self.hazard.layer.get_projection(),
            geotransform=self.hazard.layer.get_geotransform(),
            name=self.metadata().key('layer_name'),
            keywords=impact_layer_keywords,
            style_info=style_info)

        impact_layer.impact_data = impact_data
        self._impact = impact_layer
        return impact_layer

def convert_netcdf2tif(filename, n):
    """Convert netcdf to tif aggregating first n bands

    Args
        * filename: NetCDF multiband raster with extension .nc
        * n: Positive integer determining how many bands to use

    Returns
        * Raster file in tif format. Each pixel will be the maximum
          of that pixel in the first n bands in the input file.

    """

    if not isinstance(filename, basestring):
        msg = 'Argument filename should be a string. I got %s' % filename
        raise RuntimeError(msg)

    basename, ext = os.path.splitext(filename)
    msg = ('Expected NetCDF file with extension .nc - '
           'Instead I got %s' % filename)
    if ext != '.nc':
        raise RuntimeError(msg)

    try:
        n = int(n)
    except:
        msg = 'Argument N should be an integer. I got %s' % n
        raise RuntimeError(msg)

    print filename, n

    # Read NetCDF file
    fid = NetCDFFile(filename)
    dimensions = fid.dimensions.keys()
    variables = fid.variables.keys()

    title = getattr(fid, 'title')
    institution = getattr(fid, 'institution')
    source = getattr(fid, 'source')
    history = getattr(fid, 'history')
    references = getattr(fid, 'references')
    conventions = getattr(fid, 'Conventions')
    coordinate_system = getattr(fid, 'coordinate_system')

    print 'Read from %s' % filename
    print 'Title: %s' % title
    print 'Institution: %s' % institution
    print 'Source: %s' % source
    print 'History: %s' % history
    print 'References: %s' % references
    print 'Conventions: %s' % conventions
    print 'Coordinate system: %s' % coordinate_system

    print 'Dimensions: %s' % dimensions
    print 'Variables:  %s' % variables

    # Get data
    x = fid.variables['x'][:]
    y = fid.variables['y'][:]
    t = fid.variables['time'][:]
    inundation_depth = fid.variables['Inundation_Depth'][:]

    T = inundation_depth.shape[0]  # Number of time steps
    M = inundation_depth.shape[1]  # Steps in the y direction
    N = inundation_depth.shape[2]  # Steps in the x direction

    # Compute the max of the first n timesteps
    A = numpy.zeros((M, N), dtype='float')
    for i in range(n):
        B = inundation_depth[i, :, :]
        A = numpy.maximum(A, B)

    geotransform = raster_geometry2geotransform(x, y)
    print 'Geotransform', geotransform

    # Write result to tif file
    # NOTE: This assumes a default projection (WGS 84, geographic)
    date = os.path.split(basename)[-1].split('_')[0]
    print 'date', date
    R = Raster(data=A,
               geotransform=geotransform,
               keywords={'category': 'hazard',
                         'subcategory': 'flood',
                         'unit': 'm',
                         'title': ('%d hour flood forecast '
                                   'in Jakarta at %s' % (n, date))})

    tif_filename = '%s_%d_hours.tif' % (basename, n)
    R.write_to_file(tif_filename)

    print 'Success: %d hour forecast written to %s' % (n, R.filename)
    return tif_filename

#.........这里部分代码省略.........
            #    displacements > fatalities, displacements - fatalities, 0)

            # Sum up numbers for map
            # We need to use matrices here and not just numbers #2235
            # filter out NaN to avoid overflow additions
            mmi_matches = numpy.nan_to_num(mmi_matches)
            mask += mmi_matches   # Displaced

            # Generate text with result for this study
            # This is what is used in the real time system exposure table
            number_of_exposed[mmi] = exposed
            number_of_displaced[mmi] = displacements
            # noinspection PyUnresolvedReferences
            number_of_fatalities[mmi] = fatalities

        # Total statistics
        total_fatalities_raw = numpy.nansum(
            number_of_fatalities.values(), axis=0)

        # Compute probability of fatality in each magnitude bin
        if (self.__class__.__name__ == 'PAGFatalityFunction') or (
                self.__class__.__name__ == 'ITBBayesianFatalityFunction'):
            prob_fatality_mag = self.compute_probability(total_fatalities_raw)
        else:
            prob_fatality_mag = None

        # Compute number of fatalities
        self.total_population = numpy.nansum(number_of_exposed.values())
        self.total_fatalities = numpy.median(total_fatalities_raw)
        total_displaced = numpy.nansum(number_of_displaced.values())

        # As per email discussion with Ole, Trevor, Hadi, total fatalities < 50
        # will be rounded down to 0 - Tim
        # Needs to revisit but keep it alive for the time being - Hyeuk, Jono
        if self.total_fatalities < 50:
            self.total_fatalities = 0

        affected_population = self.affected_population
        affected_population[tr('Number of fatalities')] = self.total_fatalities
        affected_population[
            tr('Number of people displaced')] = total_displaced
        self.unaffected_population = (
            self.total_population - total_displaced - self.total_fatalities)
        self._evacuation_category = tr('Number of people displaced')

        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]
        total_needs = self.total_needs

        # Create style
        colours = ['#EEFFEE', '#FFFF7F', '#E15500', '#E4001B', '#730000']
        classes = create_classes(mask.flat[:], len(colours))
        interval_classes = humanize_class(classes)
        style_classes = []
        for i in xrange(len(interval_classes)):
            style_class = dict()
            style_class['label'] = create_label(interval_classes[i])
            style_class['quantity'] = classes[i]
            style_class['transparency'] = 30
            style_class['colour'] = colours[i]
            style_classes.append(style_class)

        style_info = dict(target_field=None,
                          style_classes=style_classes,
                          style_type='rasterStyle')

        impact_data = self.generate_data()

        extra_keywords = {
            'exposed_per_mmi': number_of_exposed,
            'total_population': self.total_population,
            'total_fatalities': population_rounding(self.total_fatalities),
            'total_fatalities_raw': self.total_fatalities,
            'fatalities_per_mmi': number_of_fatalities,
            'total_displaced': population_rounding(total_displaced),
            'displaced_per_mmi': number_of_displaced,
            'map_title': self.metadata().key('map_title'),
            'legend_notes': self.metadata().key('legend_notes'),
            'legend_units': self.metadata().key('legend_units'),
            'legend_title': self.metadata().key('legend_title'),
            'total_needs': total_needs,
            'prob_fatality_mag': prob_fatality_mag,
        }

        impact_layer_keywords = self.generate_impact_keywords(extra_keywords)

        # Create raster object and return
        impact_layer = Raster(
            mask,
            projection=self.exposure.layer.get_projection(),
            geotransform=self.exposure.layer.get_geotransform(),
            keywords=impact_layer_keywords,
            name=self.metadata().key('layer_name'),
            style_info=style_info)

        impact_layer.impact_data = impact_data
        self._impact = impact_layer
        return impact_layer

#.........这里部分代码省略.........
            self.volcano_names = volcano_names[:-2]  # Strip trailing ', '

        # Run interpolation function for polygon2raster
        interpolated_layer, covered_exposure_layer = \
            assign_hazard_values_to_exposure_data(
                self.hazard.layer,
                self.exposure.layer,
                attribute_name=self.target_field
            )

        # Initialise affected population per categories
        for radius in radii:
            category = 'Radius %s km ' % format_int(radius)
            self.affected_population[category] = 0

        if has_no_data(self.exposure.layer.get_data(nan=True)):
            self.no_data_warning = True
        # Count affected population per polygon and total
        for row in interpolated_layer.get_data():
            # Get population at this location
            population = row[self.target_field]
            if not numpy.isnan(population):
                population = float(population)
                # Update population count for this category
                category = 'Radius %s km ' % format_int(
                    row[self.hazard_zone_attribute])
                self.affected_population[category] += population

        # Count totals
        self.total_population = population_rounding(
            int(numpy.nansum(self.exposure.layer.get_data())))

        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]

        # Create style
        colours = ['#FFFFFF', '#38A800', '#79C900', '#CEED00',
                   '#FFCC00', '#FF6600', '#FF0000', '#7A0000']
        classes = create_classes(
            covered_exposure_layer.get_data().flat[:], len(colours))
        interval_classes = humanize_class(classes)
        # Define style info for output polygons showing population counts
        style_classes = []
        for i in xrange(len(colours)):
            style_class = dict()
            style_class['label'] = create_label(interval_classes[i])
            if i == 1:
                label = create_label(
                    interval_classes[i],
                    tr('Low Population [%i people/cell]' % classes[i]))
            elif i == 4:
                label = create_label(
                    interval_classes[i],
                    tr('Medium Population [%i people/cell]' % classes[i]))
            elif i == 7:
                label = create_label(
                    interval_classes[i],
                    tr('High Population [%i people/cell]' % classes[i]))
            else:
                label = create_label(interval_classes[i])

            style_class['label'] = label
            style_class['quantity'] = classes[i]
            style_class['colour'] = colours[i]
            style_class['transparency'] = 0
            style_classes.append(style_class)

        # Override style info with new classes and name
        style_info = dict(
            target_field=None,
            style_classes=style_classes,
            style_type='rasterStyle')

        impact_data = self.generate_data()

        # Create vector layer and return
        extra_keywords = {
            'target_field': self.target_field,
            'map_title': self.metadata().key('map_title'),
            'legend_notes': self.metadata().key('legend_notes'),
            'legend_units': self.metadata().key('legend_units'),
            'legend_title': self.metadata().key('legend_title'),
            'total_needs': self.total_needs
        }

        impact_layer_keywords = self.generate_impact_keywords(extra_keywords)

        impact_layer = Raster(
            data=covered_exposure_layer.get_data(),
            projection=covered_exposure_layer.get_projection(),
            geotransform=covered_exposure_layer.get_geotransform(),
            name=self.metadata().key('layer_name'),
            keywords=impact_layer_keywords,
            style_info=style_info)

        impact_layer.impact_data = impact_data
        self._impact = impact_layer
        return impact_layer

#.........这里部分代码省略.........
        # merely initialize
        impact = None
        for i, lo in enumerate(thresholds):
            if i == len(thresholds) - 1:
                # The last threshold
                thresholds_name = tr(
                    'People in >= %.1f m of water') % lo
                impact = medium = numpy.where(data >= lo, population, 0)
                self.impact_category_ordering.append(thresholds_name)
                self._evacuation_category = thresholds_name
            else:
                # Intermediate thresholds
                hi = thresholds[i + 1]
                thresholds_name = tr(
                    'People in %.1f m to %.1f m of water' % (lo, hi))
                medium = numpy.where((data >= lo) * (data < hi), population, 0)

            # Count
            val = int(numpy.nansum(medium))
            self.affected_population[thresholds_name] = val

        # Put the deepest area in top #2385
        self.impact_category_ordering.reverse()

        # Carry the no data values forward to the impact layer.
        impact = numpy.where(numpy.isnan(population), numpy.nan, impact)
        impact = numpy.where(numpy.isnan(data), numpy.nan, impact)

        # Count totals
        self.total_population = int(numpy.nansum(population))
        self.unaffected_population = (
            self.total_population - self.total_affected_population)

        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]

        # check for zero impact
        if numpy.nanmax(impact) == 0 == numpy.nanmin(impact):
            message = m.Message()
            message.add(self.question)
            message.add(tr('No people in %.1f m of water') % thresholds[-1])
            message = message.to_html(suppress_newlines=True)
            raise ZeroImpactException(message)

        # Create style
        colours = [
            '#FFFFFF', '#38A800', '#79C900', '#CEED00',
            '#FFCC00', '#FF6600', '#FF0000', '#7A0000']
        classes = create_classes(impact.flat[:], len(colours))
        interval_classes = humanize_class(classes)
        style_classes = []

        for i in xrange(len(colours)):
            style_class = dict()
            if i == 1:
                label = create_label(interval_classes[i], 'Low')
            elif i == 4:
                label = create_label(interval_classes[i], 'Medium')
            elif i == 7:
                label = create_label(interval_classes[i], 'High')
            else:
                label = create_label(interval_classes[i])
            style_class['label'] = label
            style_class['quantity'] = classes[i]
            style_class['transparency'] = 0
            style_class['colour'] = colours[i]
            style_classes.append(style_class)

        style_info = dict(
            target_field=None,
            style_classes=style_classes,
            style_type='rasterStyle')

        impact_data = self.generate_data()

        extra_keywords = {
            'map_title': self.metadata().key('map_title'),
            'legend_notes': self.metadata().key('legend_notes'),
            'legend_units': self.metadata().key('legend_units'),
            'legend_title': self.metadata().key('legend_title'),
            'evacuated': self.total_evacuated,
            'total_needs': self.total_needs
        }

        impact_layer_keywords = self.generate_impact_keywords(extra_keywords)

        # Create raster object and return
        impact_layer = Raster(
            impact,
            projection=self.hazard.layer.get_projection(),
            geotransform=self.hazard.layer.get_geotransform(),
            name=self.metadata().key('layer_name'),
            keywords=impact_layer_keywords,
            style_info=style_info)

        impact_layer.impact_data = impact_data
        self._impact = impact_layer
        return impact_layer

def convert_netcdf2tif(filename, n):
    """Convert netcdf to tif aggregating first n bands

    Args
        * filename: NetCDF multiband raster with extension .nc
        * n: Positive integer determining how many bands to use

    Returns
        * Raster file in tif format. Each pixel will be the maximum
          of that pixel in the first n bands in the input file.

    """

    if not isinstance(filename, basestring):
        msg = "Argument filename should be a string. I got %s" % filename
        raise RuntimeError(msg)

    basename, ext = os.path.splitext(filename)
    msg = "Expected NetCDF file with extension .nc - " "Instead I got %s" % filename
    if ext != ".nc":
        raise RuntimeError(msg)

    try:
        n = int(n)
    except:
        msg = "Argument N should be an integer. I got %s" % n
        raise RuntimeError(msg)

    print filename, n

    # Read NetCDF file
    fid = NetCDFFile(filename)
    dimensions = fid.dimensions.keys()
    variables = fid.variables.keys()

    title = getattr(fid, "title")
    institution = getattr(fid, "institution")
    source = getattr(fid, "source")
    history = getattr(fid, "history")
    references = getattr(fid, "references")
    conventions = getattr(fid, "Conventions")
    coordinate_system = getattr(fid, "coordinate_system")

    print "Read from %s" % filename
    print "Title: %s" % title
    print "Institution: %s" % institution
    print "Source: %s" % source
    print "History: %s" % history
    print "References: %s" % references
    print "Conventions: %s" % conventions
    print "Coordinate system: %s" % coordinate_system

    print "Dimensions: %s" % dimensions
    print "Variables:  %s" % variables

    # Get data
    x = fid.variables["x"][:]
    y = fid.variables["y"][:]
    t = fid.variables["time"][:]
    inundation_depth = fid.variables["Inundation_Depth"][:]

    T = inundation_depth.shape[0]  # Number of time steps
    M = inundation_depth.shape[1]  # Steps in the y direction
    N = inundation_depth.shape[2]  # Steps in the x direction

    # Compute the max of the first n timesteps
    A = numpy.zeros((M, N), dtype="float")
    for i in range(n):
        B = inundation_depth[i, :, :]
        A = numpy.maximum(A, B)

    geotransform = raster_geometry2geotransform(x, y)
    print "Geotransform", geotransform

    # Write result to tif file
    # NOTE: This assumes a default projection (WGS 84, geographic)
    date = os.path.split(basename)[-1].split("_")[0]
    print "date", date
    R = Raster(
        data=A,
        geotransform=geotransform,
        keywords={
            "category": "hazard",
            "subcategory": "flood",
            "unit": "m",
            "title": ("%d hour flood forecast grid " "in Jakarta at %s" % (n, date)),
        },
    )

    tif_filename = "%s_%d_hours.tif" % (basename, n)
    R.write_to_file(tif_filename)

    print "Success: %d hour forecast written to %s" % (n, R.filename)
    return tif_filename

#.........这里部分代码省略.........
           'Instead I got %s' % filename)
    if ext != '.nc':
        raise RuntimeError(msg)

    try:
        n = int(n)
    except:
        msg = 'Argument N should be an integer. I got %s' % n
        raise RuntimeError(msg)

    if verbose:
        print filename, n, 'hours'

    # Read NetCDF file
    fid = NetCDFFile(filename)
    dimensions = fid.dimensions.keys()
    variables = fid.variables.keys()

    title = getattr(fid, 'title')
    institution = getattr(fid, 'institution')
    source = getattr(fid, 'source')
    history = getattr(fid, 'history')
    references = getattr(fid, 'references')
    conventions = getattr(fid, 'Conventions')
    coordinate_system = getattr(fid, 'coordinate_system')

    if verbose:
        print 'Read from %s' % filename
        print 'Title: %s' % title
        print 'Institution: %s' % institution
        print 'Source: %s' % source
        print 'History: %s' % history
        print 'References: %s' % references
        print 'Conventions: %s' % conventions
        print 'Coordinate system: %s' % coordinate_system

        print 'Dimensions: %s' % dimensions
        print 'Variables:  %s' % variables

    # Get data
    x = fid.variables['x'][:]
    y = fid.variables['y'][:]
    # t = fid.variables['time'][:]
    inundation_depth = fid.variables['Inundation_Depth'][:]

    T = inundation_depth.shape[0]  # Number of time steps
    M = inundation_depth.shape[1]  # Steps in the y direction
    N = inundation_depth.shape[2]  # Steps in the x direction

    if n > T:
        msg = ('You requested %i hours prediction, but the '
               'forecast only contains %i hours' % (n, T))
        raise RuntimeError(msg)

    # Compute the max of the first n timesteps
    A = numpy.zeros((M, N), dtype='float')
    for i in range(n):
        B = inundation_depth[i, :, :]
        A = numpy.maximum(A, B)

        # Calculate overall maximal value
        total_max = numpy.max(A[:])
        #print i, numpy.max(B[:]), total_max

    geotransform = raster_geometry_to_geotransform(x, y)

    # Write result to tif file
    # NOTE: This assumes a default projection (WGS 84, geographic)
    date = os.path.split(basename)[-1].split('_')[0]

    if verbose:
        print 'Overall max depth over %i hours: %.2f m' % (n, total_max)
        print 'Geotransform', geotransform
        print 'date', date

    # Flip array upside down as it comes with rows ordered from south to north
    A = numpy.flipud(A)

    R = Raster(data=A,
               geotransform=geotransform,
               keywords={'category': 'hazard',
                         'subcategory': 'flood',
                         'unit': 'm',
                         'title': ('%d hour flood forecast grid '
                                   'in Jakarta at %s' % (n, date))})

    tif_filename = '%s_%d_hours_max_%.2f.tif' % (basename, n, total_max)
    if output_dir is not None:
        subdir_name = os.path.splitext(os.path.basename(tif_filename))[0]
        shapefile_dir = os.path.join(output_dir, subdir_name)
        if not os.path.isdir(shapefile_dir):
            os.mkdir(shapefile_dir)
        tif_filename = os.path.join(shapefile_dir, subdir_name + '.tif')

    R.write_to_file(tif_filename)

    if verbose:
        print 'Success: %d hour forecast written to %s' % (n, R.filename)

    return tif_filename

def convert_netcdf2tif(filename, n, verbose=False, output_dir=None):

    """Convert netcdf to tif aggregating first n bands

    Args
        * filename: NetCDF multiband raster with extension .nc
        * n: Positive integer determining how many bands to use
        * verbose: Boolean flag controlling whether diagnostics
          will be printed to screen. This is useful when run from
          a command line script.

    Returns
        * Raster file in tif format. Each pixel will be the maximum
          of that pixel in the first n bands in the input file.

    """

    if not isinstance(filename, basestring):
        msg = "Argument filename should be a string. I got %s" % filename
        raise RuntimeError(msg)

    basename, ext = os.path.splitext(filename)
    msg = "Expected NetCDF file with extension .nc - " "Instead I got %s" % filename
    if ext != ".nc":
        raise RuntimeError(msg)

    try:
        n = int(n)
    except:
        msg = "Argument N should be an integer. I got %s" % n
        raise RuntimeError(msg)

    if verbose:
        print filename, n, "hours"

    # Read NetCDF file
    fid = NetCDFFile(filename)
    dimensions = fid.dimensions.keys()
    variables = fid.variables.keys()

    title = getattr(fid, "title")
    institution = getattr(fid, "institution")
    source = getattr(fid, "source")
    history = getattr(fid, "history")
    references = getattr(fid, "references")
    conventions = getattr(fid, "Conventions")
    coordinate_system = getattr(fid, "coordinate_system")

    if verbose:
        print "Read from %s" % filename
        print "Title: %s" % title
        print "Institution: %s" % institution
        print "Source: %s" % source
        print "History: %s" % history
        print "References: %s" % references
        print "Conventions: %s" % conventions
        print "Coordinate system: %s" % coordinate_system

        print "Dimensions: %s" % dimensions
        print "Variables:  %s" % variables

    # Get data
    x = fid.variables["x"][:]
    y = fid.variables["y"][:]
    # t = fid.variables['time'][:]
    inundation_depth = fid.variables["Inundation_Depth"][:]

    T = inundation_depth.shape[0]  # Number of time steps
    M = inundation_depth.shape[1]  # Steps in the y direction
    N = inundation_depth.shape[2]  # Steps in the x direction

    if n > T:
        msg = "You requested %i hours prediction, but the " "forecast only contains %i hours" % (n, T)
        raise RuntimeError(msg)

    # Compute the max of the first n timesteps
    A = numpy.zeros((M, N), dtype="float")
    for i in range(n):
        B = inundation_depth[i, :, :]
        A = numpy.maximum(A, B)

        # Calculate overall maximal value
        total_max = numpy.max(A[:])
        # print i, numpy.max(B[:]), total_max

    geotransform = raster_geometry2geotransform(x, y)

    # Write result to tif file
    # NOTE: This assumes a default projection (WGS 84, geographic)
    date = os.path.split(basename)[-1].split("_")[0]

    if verbose:
        print "Overall max depth over %i hours: %.2f m" % (n, total_max)
        print "Geotransform", geotransform
        print "date", date

    # Flip array upside down as it comes with rows ordered from south to north
    A = numpy.flipud(A)

    R = Raster(
#.........这里部分代码省略.........

def convert_netcdf2tif(filename, n):
    """Convert netcdf to tif aggregating firsts n bands
    """

    if not isinstance(filename, basestring):
        msg = 'Argument filename should be a string. I got %s' % filename
        raise RuntimeError(msg)

    basename, ext = os.path.splitext(filename)
    msg = ('Expected NetCDF file with extension .nc - '
           'Instead I got %s' % filename)
    if ext != '.nc':
        raise RuntimeError(msg)

    try:
        n = int(n)
    except:
        msg = 'Argument N should be an integer. I got %s' % n
        raise RuntimeError(msg)

    print filename, n

    # Read NetCDF file
    fid = NetCDFFile(filename)
    dimensions = fid.dimensions.keys()
    variables = fid.variables.keys()

    title = getattr(fid, 'title')
    institution = getattr(fid, 'institution')
    source = getattr(fid, 'source')
    history = getattr(fid, 'history')
    references = getattr(fid, 'references')
    conventions = getattr(fid, 'Conventions')
    coordinate_system = getattr(fid, 'coordinate_system')

    print 'Read from %s' % filename
    print 'Title: %s' % title
    print 'Institution: %s' % institution
    print 'Source: %s' % source
    print 'History: %s' % history
    print 'References: %s' % references
    print 'Conventions: %s' % conventions
    print 'Coordinate system: %s' % coordinate_system

    print 'Dimensions: %s' % dimensions
    print 'Variables:  %s' % variables

    # Get data
    x = fid.variables['x'][:]
    y = fid.variables['y'][:]
    t = fid.variables['time'][:]
    inundation_depth = fid.variables['Inundation_Depth'][:]

    T = inundation_depth.shape[0]  # Number of time steps
    M = inundation_depth.shape[1]  # Steps in the y direction
    N = inundation_depth.shape[2]  # Steps in the x direction

    # Compute the max of the first n timesteps
    A = numpy.zeros((M, N), dtype='float')
    for i in range(n):
        B = inundation_depth[i, :, :]
        A = numpy.maximum(A, B)

    geotransform = raster_geometry2geotransform(x, y)
    print 'Geotransform', geotransform

    # Write result to tif file
    R = Raster(data=A,
               projection="""PROJCS["DGN95 / Indonesia TM-3 zone 48.2",
                             GEOGCS["DGN95",
                                 DATUM["Datum_Geodesi_Nasional_1995",
                                     SPHEROID["WGS 84",6378137,298.257223563,
                                         AUTHORITY["EPSG","7030"]],
                                     TOWGS84[0,0,0,0,0,0,0],
                                     AUTHORITY["EPSG","6755"]],
                                 PRIMEM["Greenwich",0,
                                     AUTHORITY["EPSG","8901"]],
                                 UNIT["degree",0.01745329251994328,
                                     AUTHORITY["EPSG","9122"]],
                                 AUTHORITY["EPSG","4755"]],
                             UNIT["metre",1,
                                 AUTHORITY["EPSG","9001"]],
                             PROJECTION["Transverse_Mercator"],
                             PARAMETER["latitude_of_origin",0],
                             PARAMETER["central_meridian",106.5],
                             PARAMETER["scale_factor",0.9999],
                             PARAMETER["false_easting",200000],
                             PARAMETER["false_northing",1500000],
                             AUTHORITY["EPSG","23834"],
                             AXIS["X",EAST],
                             AXIS["Y",NORTH]]""",
               geotransform=geotransform,
               keywords={'category': 'hazard',
                         'subcategory': 'flood',
                         'unit': 'm',
                         'title': ('Hypothetical %d hour flood forecast '
                                   'in Jakarta' % n)})
    R.write_to_file('%s_%d_hours.tif' % (basename, n))
    print 'Success: %d hour forecast written to %s' % (n, R.filename)