Python osgeo.gdal 模块，Open() 实例源码

def convert_pixgwktList_to_wgs84wktList(inputRaster, wktPolygonPixList):
    ## returns # [[GeoWKT, PixWKT], ...]
    wgs84WKTList=[]
    if os.path.isfile(inputRaster):
        srcRaster = gdal.Open(inputRaster)
        targetSR = osr.SpatialReference()
        targetSR.ImportFromWkt(srcRaster.GetProjectionRef())
        geomTransform = srcRaster.GetGeoTransform()

    for wktPolygonPix in wktPolygonPixList:
        geom_wkt_list = pixelWKTToGeoWKT(wktPolygonPix, inputRaster, targetSR='',
                                         geomTransform=geomTransform,
                                         breakMultiPolygonPix=False)

        wgs84WKTList.extend(geom_wkt_list)

    # [[GeoWKT, PixWKT], ...]
    return wgs84WKTList

def reprojectRaster(src,match_ds,dst_filename):
    src_proj = src.GetProjection()
    src_geotrans = src.GetGeoTransform()

    match_proj = match_ds.GetProjection()
    match_geotrans = match_ds.GetGeoTransform()
    wide = match_ds.RasterXSize
    high = match_ds.RasterYSize

    dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Int16)
    dst.SetGeoTransform( match_geotrans )
    dst.SetProjection( match_proj)

    gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)

    del dst # Flush
    return(gdal.Open(dst_filename,gdalconst.GA_ReadOnly))


#finds the intersecting extent of a series of scenes (left,right,bottom,top are arrays containing the respective lat/lon of the left,right,bottom,top of each image)

def get_lulc_source(ds):
    """For a given input dataset extent, select the appropriate mask source (NLCD vs. global bareground)
    """
    ds_geom = geolib.ds_geom(ds)
    lulc_fn = get_nlcd_fn()
    lulc_ds = gdal.Open(lulc_fn)
    lulc_geom = geolib.ds_geom(lulc_ds)
    #If the dem geom is within CONUS (nlcd extent), use it
    geolib.geom_transform(ds_geom, t_srs=lulc_geom.GetSpatialReference())

    if lulc_geom.Contains(ds_geom):
        print("Using NLCD 30m data for rockmask")
        lulc_source = 'nlcd'
    else:
        print("Using global 30m bare ground data for rockmask")
        #Otherwise for Global, use 30 m Bare Earth percentage 
        lulc_source = 'bareground'
        #lulc_fn = get_bareground_fn()
        #lulc_ds = gdal.Open(lulc_fn)
    return lulc_source

def read_naip(file_path, bands_to_use):
    """
    Read in a NAIP, based on www.machinalis.com/blog/python-for-geospatial-data-processing.

    Bands_to_use is an array like [0,0,0,1], designating whether to use each band (R, G, B, IR).
    """
    raster_dataset = gdal.Open(file_path, gdal.GA_ReadOnly)

    bands_data = []
    index = 0
    for b in range(1, raster_dataset.RasterCount + 1):
        band = raster_dataset.GetRasterBand(b)
        if bands_to_use[index] == 1:
            bands_data.append(band.ReadAsArray())
        index += 1
    bands_data = numpy.dstack(bands_data)

    return raster_dataset, bands_data

def tag_with_locations(test_images, predictions, tile_size, state_abbrev):
    """Combine image data with label data, so info can be rendered in a map and list UI.

    Add location data for convenience too.
    """
    combined_data = []
    for idx, img_loc_tuple in enumerate(test_images):
        raster_filename = img_loc_tuple[2]
        raster_dataset = gdal.Open(os.path.join(NAIP_DATA_DIR, raster_filename), gdal.GA_ReadOnly)
        raster_tile_x = img_loc_tuple[1][0]
        raster_tile_y = img_loc_tuple[1][1]
        ne_lat, ne_lon = pixel_to_lat_lon_web_mercator(raster_dataset, raster_tile_x +
                                                       tile_size, raster_tile_y)
        sw_lat, sw_lon = pixel_to_lat_lon_web_mercator(raster_dataset, raster_tile_x,
                                                       raster_tile_y + tile_size)
        certainty = predictions[idx][0]
        formatted_info = {'certainty': certainty, 'ne_lat': ne_lat, 'ne_lon': ne_lon,
                          'sw_lat': sw_lat, 'sw_lon': sw_lon, 'raster_tile_x': raster_tile_x,
                          'raster_tile_y': raster_tile_y, 'raster_filename': raster_filename,
                          'state_abbrev': state_abbrev, 'country_abbrev': 'USA'
                          }
        combined_data.append(formatted_info)
    return combined_data

def import_summary_geojson(geojsonfilename, removeNoBuildings=True):
    # driver = ogr.GetDriverByName('geojson')
    datasource = ogr.Open(geojsonfilename, 0)

    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())

    buildingList = []
    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:
            if removeNoBuildings:
                if feature.GetField('BuildingId') != -1:
                    buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
                                  'poly': feature.GetGeometryRef().Clone()})
            else:

                buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
                              'poly': feature.GetGeometryRef().Clone()})

    return buildingList

def import_chip_geojson(geojsonfilename, ImageId=''):
    # driver = ogr.GetDriverByName('geojson')
    datasource = ogr.Open(geojsonfilename, 0)
    if ImageId=='':
        ImageId = geojsonfilename
    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())

    polys = []
    BuildingId = 0
    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:
            BuildingId = BuildingId + 1
            polys.append({'ImageId': ImageId, 'BuildingId': BuildingId,
                          'poly': feature.GetGeometryRef().Clone()})

    return polys

def mergePolyList(geojsonfilename):

    multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
    datasource = ogr.Open(geojsonfilename, 0)

    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())


    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:

            multipolygon.AddGeometry(feature.GetGeometryRef().Clone())

    return multipolygon

def import_glcf_tile(glcf_header, cell_num, tilefile):
    glcfgrid = grids.GLCFGrid()
    glcf_data = np.zeros((glcfgrid.cell_height, glcfgrid.cell_width, 1),
                         dtype=np.uint8)

    with tempfile.NamedTemporaryFile() as f:
        # uncompress
        with gzip.open(tilefile, 'rb') as f_in, open(f.name, 'wb') as f_out:
            shutil.copyfileobj(f_in, f_out)
        gdal_ds = gdal.Open(f.name)
        assert gdal_ds is not None, "Failed to open GDAL dataset"
        band = gdal_ds.GetRasterBand(1)
        nodata_val = band.GetNoDataValue()
        print "NoData : ", nodata_val
        glcf_data[:, :, 0] = band.ReadAsArray()

    glcf_header.write_frac((cell_num, 0), glcf_data)

    print 'Finished %d' % cell_num

    return True

def reproject_tiff_from_model(old_name, new_name, model, unit):
    """
    Reprojects an tiff on a tiff model. Can be used to warp tiff.
    Keyword arguments:
    old_name -- the name of the old tiff file
    new_name -- the name of the output tiff file
    model -- the gdal dataset which will be used to warp the tiff
    unit -- the gdal unit in which the operation will be performed
    """
    mem_drv = gdal.GetDriverByName("MEM")

    old = gdal.Open(old_name)

    new = mem_drv.Create(new_name, model.RasterXSize, model.RasterYSize, 1,
                         unit)

    new.SetGeoTransform(model.GetGeoTransform())
    new.SetProjection(model.GetProjection())

    res = gdal.ReprojectImage(old, new, old.GetProjection(),
                              model.GetProjection(), gdal.GRA_NearestNeighbour)

    assert res == 0, 'Error in ReprojectImage'
    arr = new.ReadAsArray()
    new = None
    return arr

def readImageMetaData(inputFile, name):
    # Open the dataset in read only mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has been opened
    if not dataset is None:
        # Get the meta value specified
        metaData = dataset.GetMetaDataItem(name)
        # Check that it is present
        if metaData == None:
            # If not present, print error.
            print("Could not find \'", name, "\' item.")
        else:
            # Else print out the metaData value.
            print(name, "=\'", metaData, "\'")
    else:
        # Print an error messagge if the file
        # could not be opened.
        print("Could not open the input image file: ", inputFile)

# A function to read a meta-data item
# from a image band

def listImageMetaData(inputFile):
    # Open the dataset in Read Only mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has been opened.
    if not dataset is None:
        # Get the meta-data dictionary
        metaData = dataset.GetMetaData_Dict()
        # Check it has entries.
        if len(metaData) == 0:
            # if it has no entries return
            # error message.
            print("There is no image meta-data.")
        else:
            # Otherwise, print out the
            # meta-data.
            # Loop through each entry.
            for metaItem in metaData:
                print(metaItem)
    else:
        # Print an error message if the file
        # could not be opened.
        print("Could not open the input image file", inputFile)

# This is the first part of the script to
# be executed.

def setBandName(inputFile, band, name):
    # Open the image file, in update mode
    # so that the image can be edited.
    dataset = gdal.Open(inputFile, GA_Update)
    # Check that the image has opened
    if not dataset is None:
        # Get the image Band
        imgBand = dataset.GetRasterBand(band)
        # Check the image band was available
        if not imgBand is None:
            # Set the image band name
            imgBand.SetDescription(name)
        else:
            # Print out an error message
            print("Could not open the image band: ", band)

    else:
        # Print the error message if the file
        # could not be opened
        print("Could not open the input image file: ", inputFile)

# This is the first par of the script
# to be executed

def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    nodata = band.GetNoDataValue()
    array = band.ReadAsArray()

    proj = raster.GetProjection()
    inproj = osr.SpatialReference()
    inproj.ImportFromWkt(proj)

    geoTransform = raster.GetGeoTransform()
    minx = geoTransform[0]
    maxy = geoTransform[3]
    maxx = minx + geoTransform[1]*raster.RasterXSize
    miny = maxy + geoTransform[5]*raster.RasterYSize
    extent =  [minx, maxx, miny, maxy]
    pixelSizeXY = [geoTransform[1], geoTransform[5]]
    del raster, band
    return [array, nodata, extent, inproj, pixelSizeXY]

#clip a raster by vector

def gdal2np_dtype(b):
    """
    Get NumPy datatype that corresponds with GDAL RasterBand datatype
    Input can be filename, GDAL Dataset, GDAL RasterBand, or GDAL integer dtype
    """
    dt_dict = gdal_array.codes
    if isinstance(b, str):
        b = gdal.Open(b)
    if isinstance(b, gdal.Dataset):
        b = b.GetRasterBand(1)
    if isinstance(b, gdal.Band):
        b = b.DataType
    if isinstance(b, int):
        np_dtype = dt_dict[b]
    else:
        np_dtype = None
        print("Input must be GDAL Dataset or RasterBand object")
    return np_dtype

#Replace nodata value in GDAL band

def copyproj(src_fn, dst_fn, gt=True):
    """Copy projection and geotransform from one raster file to another
    """
    src_ds = gdal.Open(src_fn, gdal.GA_ReadOnly)
    dst_ds = gdal.Open(dst_fn, gdal.GA_Update)
    dst_ds.SetProjection(src_ds.GetProjection())
    if gt:
        src_gt = np.array(src_ds.GetGeoTransform())
        src_dim = np.array([src_ds.RasterXSize, src_ds.RasterYSize])
        dst_dim = np.array([dst_ds.RasterXSize, dst_ds.RasterYSize])
        #This preserves dst_fn resolution
        if np.any(src_dim != dst_dim):
            res_factor = src_dim/dst_dim.astype(float)
            src_gt[[1, 5]] *= max(res_factor)
            #src_gt[[1, 5]] *= min(res_factor)
            #src_gt[[1, 5]] *= res_factor
        dst_ds.SetGeoTransform(src_gt)
    src_ds = None
    dst_ds = None

def split_metatiles(output_dir, metatiles):

    # Determine metatile size and bit depth.
    mt_ds = gdal.Open(metatiles[0])
    assert mt_ds.RasterXSize % UF_TILE_SIZE == 0
    assert mt_ds.RasterYSize % UF_TILE_SIZE == 0

    x_uf_count = mt_ds.RasterXSize / UF_TILE_SIZE
    y_uf_count = mt_ds.RasterYSize / UF_TILE_SIZE

    for ytile in range(y_uf_count):
        print 'Process Row: ', ytile
        split_row_of_unit_fields(output_dir, metatiles, ytile)

    logging.info('Wrote %d piles to directory %s.',
                 x_uf_count * y_uf_count, output_dir)

def mask_chip(feature):
    '''
    Apply polygon mask to bounding-box chips. Chips must be named
        'feature_id.tif' and exist in the current working directory.
    '''
    fn = feature[:-4] + '.geojson'
    chip = gdal.Open(feature)

    # Create ogr vector file for gdal_rasterize
    vectordata = {'type': 'FeatureCollection', 'features': [feature]}

    with open(fn, 'wb') as f:
        geojson.dump(vectordata, f)

    # Mask raster
    cmd = 'gdal_rasterize -q -i -b 1 -b 2 -b 3 -burn 0 -burn 0 -burn 0 {} {}'.format(fn, feature)
    subprocess.call(cmd, shell=True)

    # Remove ogr vector file
    os.remove(fn)

def rasterize(self,inRaster,inShape,inField):
        filename = tempfile.mktemp('.tif')        
        data = gdal.Open(inRaster,gdal.GA_ReadOnly)
        shp = ogr.Open(inShape)

        lyr = shp.GetLayer()

        driver = gdal.GetDriverByName('GTiff')
        dst_ds = driver.Create(filename,data.RasterXSize,data.RasterYSize, 1,gdal.GDT_Byte)
        dst_ds.SetGeoTransform(data.GetGeoTransform())
        dst_ds.SetProjection(data.GetProjection())
        OPTIONS = 'ATTRIBUTE='+inField
        gdal.RasterizeLayer(dst_ds, [1], lyr, None,options=[OPTIONS])
        data,dst_ds,shp,lyr=None,None,None,None


        return filename

def clip(self):

        """
        This function clips the mosaicked, projected images to the size of the
        referenceImage

        """

        subprocess.call(['gdaltindex', self.extent, self.referenceImagePath])
        dataNames = sorted(glob.glob(self.fullPath + '/full*.tif'))
        splitAt = len(self.fullPath) + 1

        for i in range(len(dataNames)):
            x = dataNames[i]
            y = dataNames[i][:splitAt] + dataNames[i][splitAt+4:]
            subprocess.call(['gdalwarp', '-r', 'near', '-cutline', self.extent, '-crop_to_cutline', x, y, '-dstnodata', '9999'])

        for n in dataNames:
            os.remove(n)
        dataNames = sorted(glob.glob(self.fullPath + '/*.tif'))
        test = gdal.Open(dataNames[0]).ReadAsArray()
        logger.log('SUCCESS', 'Clipping complete!  %d %s files  were successfully clipped to the size of %s with dimensions %d rows by %d columns' % (len(dataNames), str(self.outformat), str(self.referenceImagePath), test.shape[0], test.shape[1]))

def rasterize_polygons(polygon_shp, template_raster, out_raster, field):
    """ generate a categorical raster based on polygons

    :rtype: None
    :param polygon_shp: input polygon shapefile
    :param template_raster: raster template for cellsize and extent
    :param out_raster: output raster file
    """

    # Source: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

    gdal.UseExceptions()
    # Get template raster specs
    src_ds = gdal.Open(template_raster)
    if src_ds is None:
        print 'Unable to open %s' % template_raster
        sys.exit(1)
    try:
        srcband = src_ds.GetRasterBand(1)
    except RuntimeError, e:
        print 'No band %i found' % 0
        print e
        sys.exit(1)

    # Open the data source and read in the extent
    source_ds = ogr.Open(polygon_shp)
    source_layer = source_ds.GetLayer()

    target_ds = gdal.GetDriverByName('GTiff').Create(out_raster, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Float32)
    target_ds.SetGeoTransform(src_ds.GetGeoTransform())
    target_ds.SetProjection(src_ds.GetProjection())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(srcband.GetNoDataValue())

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE={}".format(field)])

def readTodayBands(path,row):
    allFiles=sorted(glob.glob(os.path.join(getCurrentDirectory(),'Today',path,row,'*.TIF')))
    allRasters=dict([])
    fileNumber=1
    percent=(fileNumber/len(allFiles))*100
    for file in allFiles:
        percent=(fileNumber/len(allFiles))*100
        sys.stdout.write("\rReading today's bands...%d%%" % percent)
        sys.stdout.flush()
        bandName=file[file.rfind('_')+1:-4]
        sys.stdout.write("\rReading today's bands...%d%%" % percent)
        sys.stdout.flush()
        allRasters[bandName]=gdal.Open(file,gdalconst.GA_ReadOnly)
        fileNumber+=1
    todayExtent = getRasterExtent(allRasters['B4'])#saves the extent of today's rasters (they'll match band 4) so that we can crop during the cloudmask backfill
    print('')
    return(allRasters,todayExtent)

#large function that will back-fill cloudy areas of most recent imagery using historic imagery

def downloadLandsatScene(jd,year,month,day,path,row,directory):
    LandsatID=LSUniqueID(path,row,year,jd,os.path.join(directory,str(row)))
    downLoadCommand=getCurrentDirectory() + "/download_landsat_scene.py -o scene -b LC8 -d "+ year + month + day + ' -s ' + path + '0'+str(row)+ " -u usgs.txt --output " + LandsatID
    os.system(downLoadCommand)
    print('extracting...')
    tarName=extractTar(LandsatID,os.path.join(directory,str(row)))
    allFiles=glob.glob(os.path.join(directory,str(row),'*.TIF'))
    for filename in allFiles:
        bandName=filename.strip('.TIF')[filename.rfind('_')+1:]
        if bandName != 'B4' and bandName != 'B5' and bandName != 'B7' and bandName != 'B8' and bandName != 'BQA':
            os.remove(filename)
    try:
        shutil.rmtree(os.path.join(directory,str(row),tarName))
    except:
        print('No folder to delete called: ' + os.path.join(directory,str(row),tarName))
    try:
        os.remove(os.path.join(directory,str(row),tarName + '_MTL.txt'))
    except:
        print('No file to delete called: ' + os.path.join(directory,str(row),tarName + '_MTL.txt'))
    reprojectPanBand(gdal.Open(os.path.join(directory,str(row),tarName + '_B8.TIF'),gdalconst.GA_ReadOnly),gdal.Open(os.path.join(directory,str(row),tarName + '_B4.TIF'),gdalconst.GA_ReadOnly),os.path.join(directory,str(row),tarName + '_B8.TIF'))
    SLIP.model(year+month+day,path,str(row))

def test_hall_rectification(self):
        '''Should rectify an image in the expected way.'''
        control_sets = {
            'High/Bright': [(331501.45,4694346.66), (319495.39,4706820.66), (298527.006,4691417.99)],
            'Low/Dark': [(322577.40,4658508.99), (361612.79,4694665.62), (378823.69,4692132.56)]
        }
        ref_raster = gdal.Open(os.path.join(self.test_dir, 'multi7_raster.tiff'))
        sub_raster = gdal.Open(os.path.join(self.test_dir, 'multi7_raster2.tiff'))

        # NOTE: Using same control sets for reference, subject
        hall_rectification(ref_raster, sub_raster, self.test_dir, control_sets, control_sets)

        arr, gt, wkt = as_array(os.path.join(self.test_dir, 'rect_multi7_raster2.tiff'))
        self.assertTrue(np.array_equal(arr.shape, (6, 74, 81)))
        self.assertTrue(np.array_equiv(arr[:,50,50].round(5), np.array([
            17, 1331, 1442, 3422, 2916, 2708
        ]).round(5)))

def test_spectral_profile(self):
        '''
        Should correctly retrieve a spectral profile from a raster dataset.
        '''
        coords = ((-84.5983, 42.7256), (-85.0807, 41.1138))
        pixels = [(18, 0), (2, 59)]
        file_path = os.path.join(self.test_dir, 'multi3_raster.tiff')
        ds = gdal.Open(file_path)
        kwargs = {
            'gt': ds.GetGeoTransform(),
            'wkt': ds.GetProjection(),
            'dd': True
        }

        # The true spectral profile
        spectra = np.array([[237, 418, 325], [507, 616, 445]], dtype=np.int16)
        sp1 = spectra_at_xy(ds, coords, **kwargs)
        sp2 = spectra_at_xy(ds.ReadAsArray(), coords, **kwargs)
        sp3 = spectra_at_idx(ds.ReadAsArray().transpose(), pixels)
        self.assertEqual(spectra.tolist(), sp1.tolist())
        self.assertEqual(spectra.tolist(), sp2.tolist())
        self.assertEqual(spectra.tolist(), sp3.tolist())

def array_to_raster_clone(a, proto, xoff=None, yoff=None):
    '''
    Creates a raster from a given array and a prototype raster dataset, with
    optional x- and y-offsets if the array was clipped. Arguments:
        a       A NumPy array
        proto   A prototype dataset
        xoff    The offset in the x-direction; should be provided when clipped
        yoff    The offset in the y-direction; should be provided when clipped
    '''
    rast = gdal_array.OpenNumPyArray(a)

    kwargs = dict()
    if xoff is not None and yoff is not None:
        kwargs = dict(xoff=xoff, yoff=yoff)

    # Copy the projection info and metadata from a prototype dataset
    if type(proto) == str:
        proto = gdal.Open(proto)

    gdalnumeric.CopyDatasetInfo(proto, rast, **kwargs)

    return rast

def get_lulc_ds_full(ds, lulc_source=None):
    if lulc_source is None:
        lulc_source = get_lulc_source(ds)
    if lulc_source == 'nlcd':
        lulc_ds_full = gdal.Open(get_nlcd_fn())
    elif lulc_source == 'bareground':
        lulc_ds_full = gdal.Open(get_bareground_fn())
    return lulc_ds_full

def proc_modscag(fn_list, extent=None, t_srs=None):
    """Process the MODSCAG products for full date range, create composites and reproject
    """
    #Use cubic spline here for improve upsampling 
    ds_list = warplib.memwarp_multi_fn(fn_list, res='min', extent=extent, t_srs=t_srs, r='cubicspline')
    stack_fn = os.path.splitext(fn_list[0])[0] + '_' + os.path.splitext(os.path.split(fn_list[-1])[1])[0] + '_stack_%i' % len(fn_list) 
    #Create stack here - no need for most of mastack machinery, just make 3D array
    #Mask values greater than 100% (clouds, bad pixels, etc)
    ma_stack = np.ma.array([np.ma.masked_greater(iolib.ds_getma(ds), 100) for ds in np.array(ds_list)], dtype=np.uint8)

    stack_count = np.ma.masked_equal(ma_stack.count(axis=0), 0).astype(np.uint8)
    stack_count.set_fill_value(0)
    stack_min = ma_stack.min(axis=0).astype(np.uint8)
    stack_min.set_fill_value(0)
    stack_max = ma_stack.max(axis=0).astype(np.uint8)
    stack_max.set_fill_value(0)
    stack_med = np.ma.median(ma_stack, axis=0).astype(np.uint8)
    stack_med.set_fill_value(0)

    out_fn = stack_fn + '_count.tif'
    iolib.writeGTiff(stack_count, out_fn, ds_list[0])
    out_fn = stack_fn + '_max.tif'
    iolib.writeGTiff(stack_max, out_fn, ds_list[0])
    out_fn = stack_fn + '_min.tif'
    iolib.writeGTiff(stack_min, out_fn, ds_list[0])
    out_fn = stack_fn + '_med.tif'
    iolib.writeGTiff(stack_med, out_fn, ds_list[0])

    ds = gdal.Open(out_fn)
    return ds

def _open(self):
            if not _gdal:
                load_lib()
            self._ds = _gdal.Open(self.request.get_local_filename())

def polygonize(self,rasterTemp,outShp):

            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)

            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile

            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)

            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  

            outDatasource.Destroy()
            sourceRaster=None
            band=None


            ioShpFile = ogr.Open(outShp, update = 1)


            lyr = ioShpFile.GetLayerByIndex(0)

            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()

            return outShp

def test_gdal_vs_numpy(self):
        for t in self.tifs:
            ds = gdal.Open(t)
            n_list = geoinfo.numpy_band_stats(ds, t, 1)
            g_list = geoinfo.band_stats(ds, t, 1)
            self.assertEqual(n_list[0], g_list[0])
            self.assertEqual(n_list[-3], g_list[-3])
            self.assertEqual(n_list[-2], g_list[-2])

            np.testing.assert_almost_equal(
                np.array([float(t) for t in n_list[1:-3]]),
                np.array([float(t) for t in g_list[1:-3]]),
                decimal=4
            )

def test_geotransform_projection_nodata(self):
        tmp_output = tempfile.mktemp(suffix='.tif')
        extents = [str(s) for s in [920000, -4300000, 929000, -4290000]]

        # the mock is already geotransformed, so this will have no effect
        # to projection and nodata, but will be cropped making the
        # geotransform tuple different
        crop.crop_reproject_resample(self.std2000_no_mask, tmp_output,
                                     sampling='bilinear',
                                     extents=extents,
                                     reproject=True)

        ds = gdal.Open(tmp_output)
        gt = ds.GetGeoTransform()
        projection = ds.GetProjection()
        nodata = ds.GetRasterBand(1).GetNoDataValue()
        ds = None

        ds = gdal.Open(self.std2000_no_mask)
        projection_input = ds.GetProjection()
        nodata_input = ds.GetRasterBand(1).GetNoDataValue()
        ds = None

        self.assertEqual(nodata, nodata_input)
        self.assertEqual(projection, projection_input)
        self.assertEqual(gt[1], float(crop.OUTPUT_RES[0]))
        self.assertEqual(gt[0], float(extents[0]))
        self.assertEqual(gt[3], float(extents[3]))
        os.remove(tmp_output)

def test_apply_mask(self):
        output_file = tempfile.mktemp(suffix='.tif')
        jpeg = False
        tmp_out_file = tempfile.mktemp(suffix='.tif')
        shutil.copy(self.std2000_no_mask, tmp_out_file)
        crop.apply_mask(mask_file=self.mask,
                        tmp_output_file=tmp_out_file,
                        output_file=output_file,
                        jpeg=jpeg)

        ds = gdal.Open(output_file)
        gt = ds.GetGeoTransform()
        projection = ds.GetProjection()
        nodata = ds.GetRasterBand(1).GetNoDataValue()
        ds = None

        ds = gdal.Open(self.std2000)
        projection_input = ds.GetProjection()
        nodata_input = ds.GetRasterBand(1).GetNoDataValue()
        ds = None

        self.assertEqual(nodata, nodata_input)
        self.assertEqual(projection, projection_input)
        self.assertEqual(gt[1], float(crop.OUTPUT_RES[0]))
        self.assertEqual(gt[0], self.extents[0])
        self.assertEqual(gt[3], self.extents[3])
        os.remove(output_file)

def test_mpi_vs_serial(self):

        tmpdir1 = tempfile.mkdtemp()
        tmpdir2 = tempfile.mkdtemp()

        for n, partitions in product(range(1, 3), range(2, 100, 10)):
            size = 2 * n + 1
            raster_average.treat_file(self.test_tif,
                                      out_dir=tmpdir1,
                                      size=size,
                                      func='nanmean',
                                      partitions=partitions)
            arr1 = gdal.Open(
                join(tmpdir1, basename(self.test_tif))).ReadAsArray()

            raster_average.treat_file(self.test_tif,
                                      out_dir=tmpdir2,
                                      size=size,
                                      func='nanmean',
                                      partitions=partitions)
            arr2 = gdal.Open(join(tmpdir2,
                                  basename(self.test_tif))).ReadAsArray()

            np.testing.assert_array_almost_equal(arr1, arr2)

        shutil.rmtree(tmpdir2)
        shutil.rmtree(tmpdir1)

def apply_mask(mask_file, tmp_output_file, output_file, jpeg):
    """
    Parameters
    ----------
    mask_file: mask file path
    tmp_output_file: intermediate cropped geotiff before mask application
    output_file: output geotiff path
    jpeg: boolean, whether to produce jpeg or not
    -------

    """
    mask = get_mask(mask_file)

    out_ds = gdal.Open(tmp_output_file, gdal.GA_Update)
    out_band = out_ds.GetRasterBand(1)
    out_data = out_band.ReadAsArray()
    no_data_value = out_band.GetNoDataValue()
    log.info('Found NoDataValue {} for file {}'.format(
        no_data_value, os.path.basename(tmp_output_file)))
    if no_data_value is not None:
        out_data[mask] = no_data_value
    else:
        log.warning('NoDataValue was not set for {}'.format(tmp_output_file))
        log.info('Manually setting NoDataValue for {}'.format(tmp_output_file))
    out_band.WriteArray(out_data)
    out_ds = None  # close dataset and flush cache

    # move file to output file
    shutil.move(tmp_output_file, output_file)
    log.info('Output file {} created'.format(tmp_output_file))

    if jpeg:
        dir_name = os.path.dirname(output_file)
        jpeg_file = os.path.basename(output_file).split('.')[0] + '.jpg'
        jpeg_file = os.path.join(dir_name, jpeg_file)
        cmd_jpg = ['gdal_translate', '-ot', 'Byte', '-of', 'JPEG', '-scale',
                   output_file,
                   jpeg_file] + COMMON
        subprocess.check_call(cmd_jpg)
        log.info('Created {}'.format(jpeg_file))

def average(input_dir, out_dir, size):
    input_dir = abspath(input_dir)
    log.info('Reading tifs from {}'.format(input_dir))
    tifs = glob.glob(join(input_dir, '*.tif'))

    for tif in tifs:
        data_source = gdal.Open(tif, gdal.GA_ReadOnly)
        band = data_source.GetRasterBand(1)
        # data_type = gdal.GetDataTypeName(band.DataType)
        data = band.ReadAsArray()
        no_data_val = band.GetNoDataValue()
        averaged_data = filter_data(data, size, no_data_val)
        log.info('Calculated average for {}'.format(basename(tif)))

        output_file = join(out_dir, 'average_' + basename(tif))
        out_ds = gdal.GetDriverByName('GTiff').Create(
            output_file, data_source.RasterXSize, data_source.RasterYSize,
            1, band.DataType)
        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(averaged_data)
        out_ds.SetGeoTransform(data_source.GetGeoTransform())
        out_ds.SetProjection(data_source.GetProjection())
        out_band.FlushCache()  # Write data to disc
        out_ds = None  # close out_ds
        data_source = None  # close dataset

        log.info('Finished converting {}'.format(basename(tif)))

def gdalaverage(input_dir, out_dir, size):
    """
    average data using gdal's averaging method.
    Parameters
    ----------
    input_dir: str
        input dir name of the tifs that needs to be averaged
    out_dir: str
        output dir name
    size: int, optional
        size of kernel
    Returns
    -------

    """
    input_dir = abspath(input_dir)
    log.info('Reading tifs from {}'.format(input_dir))
    tifs = glob.glob(join(input_dir, '*.tif'))

    process_tifs = np.array_split(tifs, mpiops.chunks)[mpiops.chunk_index]

    for tif in process_tifs:
        data_set = gdal.Open(tif, gdal.GA_ReadOnly)
        # band = data_set.GetRasterBand(1)
        # data_type = gdal.GetDataTypeName(band.DataType)
        # data = band.ReadAsArray()
        # no_data_val = band.GetNoDataValue()
        # averaged_data = filter_data(data, size, no_data_val)
        log.info('Calculated average for {}'.format(basename(tif)))

        output_file = join(out_dir, 'average_' + basename(tif))
        src_gt = data_set.GetGeoTransform()
        tmp_file = '/tmp/tmp_{}.tif'.format(mpiops.chunk_index)
        resample_cmd = [TRANSLATE] + [tif, tmp_file] + \
            ['-tr', str(src_gt[1]*size), str(src_gt[1]*size)] + \
            ['-r', 'bilinear']
        check_call(resample_cmd)
        rollback_cmd = [TRANSLATE] + [tmp_file, output_file] + \
            ['-tr', str(src_gt[1]), str(src_gt[1])]
        check_call(rollback_cmd)
        log.info('Finished converting {}'.format(basename(tif)))

def get_stats(tif, partitions=100):
    ds = gdal.Open(tif, gdal.GA_ReadOnly)
    number_of_bands = ds.RasterCount

    if ds.RasterCount > 1:
        d = {}
        log.info('Found multibanded geotif {}'.format(basename(tif)))
        for b in range(number_of_bands):
            d['{tif}_{b}'.format(tif=tif, b=b)] = band_stats(ds, tif, b + 1,
                                                             partitions)
        return d
    else:
        log.info('Found single band geotif {}'.format(basename(tif)))
        return {tif: band_stats(ds, tif, 1, partitions)}

def get_numpy_stats(tif, writer):
    ds = gdal.Open(tif, gdal.GA_ReadOnly)
    number_of_bands = ds.RasterCount

    if ds.RasterCount > 1:
        log.info('Found multibanded geotif {}'.format(basename(tif)))
        for b in range(number_of_bands):
            write_rows(stats=number_of_bands(ds, tif, b + 1), writer=writer)
    else:
        log.info('Found single band geotif {}'.format(basename(tif)))
        write_rows(stats=number_of_bands(ds, tif, 1), writer=writer)

def VPMprmt (directory,tile):
    tempfile=None
    for path, subdirs, files in os.walk(directory):
        for name in files:
            if (tile in name) and (".tif" == name[-4:]) : tempfile=os.path.join(path,name)
    print tempfile
    inprmtdata=gdal.Open(tempfile)
    prmtdata=inprmtdata.ReadAsArray()
    return prmtdata


#Function to build vrt

def movingaverage(tile):
    # Output directories for moving average LSWImax

    # if the output directories don't exist, create the new directories
    if not os.path.exists(dirLSWIMA):
        os.makedirs(dirLSWIMA)
    # build LSWImax vrt file and read as an array
    file=buildVrtFile(year,tile)
    vrtLSWImax=os.path.join(os.path.dirname(file),year+tile+'LSWImax_vrt.vrt')
    #print "Building the vrt file: ", year+tile+vrtLSWImax
    os.system('gdalbuildvrt -separate -input_file_list '+file+' '+vrtLSWImax)

    global rows, cols, geoProj,geoTran
    inLSWImax=gdal.Open(vrtLSWImax)
    #print "reading the multi-LSWI..."
    LSWImax=inLSWImax.ReadAsArray()

    rows = 2400
    cols = 2400
    geoTran=inLSWImax.GetGeoTransform()
    geoProj=inLSWImax.GetProjection()

    #find the second largest LSWImax
    secLSWImax = numpy.sort(LSWImax, axis=0, kind='quicksort', order=None)[3,:,:]

    write_file(dirLSWIMA+'/'+tile+'.'+year+'.maxLSWI_MA.tif',secLSWImax,geoTran,rows,cols,geoProj,driverName='GTiff')

    secLSWImax=None
    LSWImax=None

def movingaverage(tile):
    # Output directories for moving average LSWImax

    # if the output directories don't exist, create the new directories
    if not os.path.exists(dirLSWIMA):
        os.makedirs(dirLSWIMA)
    # build LSWImax vrt file and read as an array
    file=buildVrtFile(year,tile)
    vrtLSWImax=os.path.join(os.path.dirname(file),year+tile+'LSWImax_vrt.vrt')
    #print "Building the vrt file: ", year+tile+vrtLSWImax
    os.system('gdalbuildvrt -separate -input_file_list '+file+' '+vrtLSWImax)

    global rows, cols, geoProj,geoTran
    inLSWImax=gdal.Open(vrtLSWImax)
    #print "reading the multi-LSWI..."
    LSWImax=inLSWImax.ReadAsArray()

    rows = 2400
    cols = 2400
    geoTran=inLSWImax.GetGeoTransform()
    geoProj=inLSWImax.GetProjection()

    #find the second largest LSWImax
    secLSWImax = numpy.sort(LSWImax, axis=0, kind='quicksort', order=None)[3,:,:]

    write_file(dirLSWIMA+'/'+tile+'.'+year+'.maxLSWI_MA.tif',secLSWImax,geoTran,rows,cols,geoProj,driverName='GTiff')

    secLSWImax=None
    LSWImax=None

def buildTindex(rasterFolder, rasterExtention='.tif'):
    rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
    print(rasterList)

    print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))

    memDriver = ogr.GetDriverByName('MEMORY')
    gTindex = memDriver.CreateDataSource('gTindex')
    srcImage = gdal.Open(rasterList[0])
    spat_ref = osr.SpatialReference()
    spat_ref.SetProjection(srcImage.GetProjection())
    gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("location", ogr.OFTString)
    gTindexLayer.CreateField(idField)

    # Create the feature and set values
    featureDefn = gTindexLayer.GetLayerDefn()



    for rasterFile in rasterList:
        srcImage = gdal.Open(rasterFile)

        geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)

        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(polyToCut)
        feature.SetField("location", rasterFile)
        gTindexLayer.CreateFeature(feature)
        feature = None


    return gTindex, gTindexLayer

def createTiledGeoJsonFromSrc(rasterFolderLocation, vectorSrcFile, geoJsonOutputDirectory, rasterTileIndex='',
                              geoJsonPrefix='GEO', rasterFileExtenstion='.tif',
                              rasterPrefixToReplace='PAN'
                              ):
    if rasterTileIndex == '':
        gTindex, gTindexLayer = buildTindex(rasterFolderLocation, rasterExtention=rasterFileExtenstion)
    else:
        gTindex = ogr.Open(rasterTileIndex,0)
        gTindexLayer = gTindex.GetLayer()

    shapeSrc = ogr.Open(vectorSrcFile,0)
    chipSummaryList = []
    for feature in gTindexLayer:
        featureGeom = feature.GetGeometryRef()
        rasterFileName = feature.GetField('location')
        rasterFileBaseName = os.path.basename(rasterFileName)
        outGeoJson = rasterFileBaseName.replace(rasterPrefixToReplace, geoJsonPrefix)
        outGeoJson = outGeoJson.replace(rasterFileExtenstion, '.geojson')
        outGeoJson = os.path.join(geoJsonOutputDirectory, outGeoJson)

        gT.clipShapeFile(shapeSrc, outGeoJson, featureGeom, minpartialPerc=0.0, debug=False)
        imageId = rasterFileBaseName.replace(rasterPrefixToReplace+"_", "")
        chipSummary = {'chipName': rasterFileName,
                           'geoVectorName': outGeoJson,
                           'imageId': os.path.splitext(imageId)[0]}

        chipSummaryList.append(chipSummary)

    return chipSummaryList

def createTruthPixelLinePickle(truthLineFile, pickleLocation=''):
    if pickleLocation=='':
        extension = os.path.splitext(truthLineFile)[1]
        pickleLocation = truthLineFile.replace(extension, 'Pixline.p')
    if truthLineFile != '':
        # get Source Line File Information
        shapef = ogr.Open(truthLineFile, 0)
        truthLayer = shapef.GetLayer()
        pt1X = []
        pt1Y = []
        pt2X = []
        pt2Y = []
        for tmpFeature in truthLayer:
            tmpGeom = tmpFeature.GetGeometryRef()
            for i in range(0, tmpGeom.GetPointCount()):
                pt = tmpGeom.GetPoint(i)

                if i == 0:
                    pt1X.append(pt[0])
                    pt1Y.append(pt[1])
                elif i == 1:
                    pt2X.append(pt[0])
                    pt2Y.append(pt[1])

        lineData = {'pt1X': np.asarray(pt1X),
                    'pt1Y': np.asarray(pt1Y),
                    'pt2X': np.asarray(pt2X),
                    'pt2Y': np.asarray(pt2Y)
                    }

        with open(pickleLocation, 'wb') as f:
            pickle.dump(lineData, f)
            # get Source Line File Information

def createTruthPixelPolyPickle(truthPoly, pickleLocation=''):
    # returns dictionary with list of minX, maxX, minY, maxY

    if pickleLocation=='':
        extension = os.path.splitext(truthPoly)[1]
        pickleLocation = truthPoly.replace(extension, 'PixPoly.p')
    if truthPoly != '':
        # get Source Line File Information
        shapef = ogr.Open(truthPoly, 0)
        truthLayer = shapef.GetLayer()
        envList = []

        for tmpFeature in truthLayer:
            tmpGeom = tmpFeature.GetGeometryRef()
            env = tmpGeom.GetEvnelope()
            envList.append(env)

        envArray = np.asarray(envList)
        envelopeData = {'minX': envArray[:,0],
                        'maxX': envArray[:,1],
                        'minY': envArray[:,2],
                        'maxY': envArray[:,3]
                        }


        with open(pickleLocation, 'wb') as f:
            pickle.dump(envelopeData, f)
            # get Source Line File Information

def openm3(input_data):
    if input_data.split('.')[-1] == 'hdr':
        # GDAL wants the img, but many users aim at the .hdr
        input_data = input_data.split('.')[0] + '.img'
    ds = gdal.Open(input_data)
    ref_array = ds.GetRasterBand(1).ReadAsArray()
    metadata = ds.GetMetadata()
    wv_array = metadatatoband(metadata)
    return wv_array, ref_array, ds

def openmi(input_data):
    ds = gdal.Open(input_data)
    band_pointers = []
    nbands = ds.RasterCount

    for b in xrange(1, nbands + 1):
        band_pointers.append(ds.GetRasterBand(b))

    ref_array = ds.GetRasterBand(1).ReadAsArray()
    wv_array = None
    return wv_array, ref_array[::3, ::3], ds

def readBandMetaData(inputFile, band, name):
    # Open the dataset in Read Only mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has opened
    if not dataset is None:
        # Get the image band
        imgBand = dataset.GetRasterBand(band)
        # Check the image band was available
        if not imgBand is None:
            # Get the meta data value specified
            metaData = imgBand.GetMetaDataItem(name)
            # Check that it is present
            if metaData == None:
                # If not present, print error.
                print("Could not find \'", name, "\' item.")
            else:
                # Else print out the metadata value
                print(name, "=\'", metaData, "\'")

        else:
            # Print out an error message.
            print("Could not open the image band: ", band)

    else:
        # Print an error message if the file
        # could not be opened.
        print("Could not open the input image file: ", inputFile)

# A function to read a meta-data item
# from an image