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

def gdt2pt(gdt):
    """Translate GDAL data type to WKT Raster pixel type."""
    pixtypes = {
        gdalc.GDT_Byte    : { 'name': 'PT_8BUI',  'id':  4 },
        gdalc.GDT_Int16   : { 'name': 'PT_16BSI', 'id':  5 },
        gdalc.GDT_UInt16  : { 'name': 'PT_16BUI', 'id':  6 },
        gdalc.GDT_Int32   : { 'name': 'PT_32BSI', 'id':  7 },
        gdalc.GDT_UInt32  : { 'name': 'PT_32BUI', 'id':  8 },
        gdalc.GDT_Float32 : { 'name': 'PT_32BF',  'id': 10 },
        gdalc.GDT_Float64 : { 'name': 'PT_64BF',  'id': 11 }
        }

    # XXX: Uncomment these logs to debug types translation
    #logit('MSG: Input GDAL pixel type: %s (%d)\n' % (gdal.GetDataTypeName(gdt), gdt))
    #logit('MSG: Output WKTRaster pixel type: %(name)s (%(id)d)\n' % (pixtypes.get(gdt, 13)))

    return pixtypes.get(gdt, 13)

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 rasterToJSON (infile, outfile, outproj):
    #reproject dataset to the desired projection
    subprocess.call(['gdalwarp', '-t_srs', 'EPSG:{}'.format(outproj), infile, 'outReproj.tif', '-overwrite'])
    projRas = gdal.Open('outReproj.tif')

    #get properties of new raster
    nrows = projRas.RasterYSize
    ncols = projRas.RasterXSize
    band1 = projRas.GetRasterBand(1)
    nodata = band1.GetNoDataValue()
    gdata = band1.ReadAsArray()
    x0, y0 = projRas.GetGeoTransform()[0], projRas.GetGeoTransform()[3]
    cellX, cellY = projRas.GetGeoTransform()[1], projRas.GetGeoTransform()[5]

    #get corners
    ulcorner, llcorner = [x0, y0], [x0, y0 + nrows*cellY]
    urcorner, lrcorner =  [x0 + ncols*cellX, y0], [x0 + ncols*cellX, y0 + nrows*cellY]

    dataType = gdal.GetDataTypeName(band1.DataType)
    #convert to float
    if dataType.startswith('Int'):
        gdata = gdata.astype(numpy.float32, copy=False)
    # new nodata value
    gdata[gdata == nodata] = -9999
    gdata = numpy.concatenate(gdata)
    gdata = gdata.tolist()
    #write to json
    gdataJSON = {'upLeft': transPoint(ulcorner, outproj, 4326), 'loLeft': transPoint(llcorner, outproj, 4326),
                 'upRight': transPoint(urcorner, outproj, 4326), 'loRight': transPoint(lrcorner, outproj, 4326),
                 'projEPSG': outproj, 'width': ncols, 'height': nrows, 'data': gdata}

    with open(outfile, 'w') as fp:
        json.dump(gdataJSON, fp)

    del gdata
    subprocess.call(['rm', 'outReproj.tif'])

def changeRasterValues(self, band):

        fmttypes = {'Byte': 'B', 'UInt16': 'H', 'Int16': 'h', 'UInt32': 'I', 'Int32': 'i', 'Float32': 'f',
                    'Float64': 'd'}

        data_type = band.DataType

        BandType = gdal.GetDataTypeName(band.DataType)

        raster = []

        for y in range(band.YSize):
            scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, data_type)
            values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
            raster.append(values)

        raster = [list(item) for item in raster]

        # changing raster values
        #for i, item in enumerate(raster):
        #    for j, value in enumerate(item):
        #        #if value > 98:
        #        #    #print i, j
        #        #    raster[i][j] = 0

        # transforming list in array
        raster = numpy.asarray(raster)

        return raster

def __init__(self, filename):
        if os.path.isfile(filename):
            self.filename = filename if os.path.isabs(filename) else os.path.join(os.getcwd(), filename)
            self.raster = gdal.Open(filename, GA_ReadOnly)
        else:
            raise IOError('file does not exist')

        self.cols = self.raster.RasterXSize
        self.rows = self.raster.RasterYSize
        self.bands = self.raster.RasterCount
        self.dim = [self.rows, self.cols, self.bands]
        self.driver = self.raster.GetDriver()
        self.format = self.driver.ShortName
        self.dtype = gdal.GetDataTypeName(self.raster.GetRasterBand(1).DataType)
        self.projection = self.raster.GetProjection()
        self.srs = osr.SpatialReference(wkt=self.projection)
        self.proj4 = self.srs.ExportToProj4()
        try:
            self.epsg = crsConvert(self.srs, 'epsg')
        except RuntimeError:
            self.epsg = None
        self.geogcs = self.srs.GetAttrValue('geogcs')
        self.projcs = self.srs.GetAttrValue('projcs') if self.srs.IsProjected() else None
        self.geo = dict(
            zip(['xmin', 'xres', 'rotation_x', 'ymax', 'rotation_y', 'yres'], self.raster.GetGeoTransform()))

        # note: yres is negative!
        self.geo['xmax'] = self.geo['xmin'] + self.geo['xres'] * self.cols
        self.geo['ymin'] = self.geo['ymax'] + self.geo['yres'] * self.rows

        self.res = [abs(float(self.geo['xres'])), abs(float(self.geo['yres']))]
        self.nodata = self.raster.GetRasterBand(1).GetNoDataValue()

        self.__data = []

def get_ndv_b(b):
    """Get NoData value for GDAL band.

    If NoDataValue is not set in the band, 
    extract upper left and lower right pixel values.
    Otherwise assume NoDataValue is 0.

    Parameters
    ----------
    b : GDALRasterBand object 
        This is the input band.

    Returns
    -------
    b_ndv : float 
        NoData value 
    """

    b_ndv = b.GetNoDataValue()
    if b_ndv is None:
        #Check ul pixel for ndv
        ns = b.XSize
        nl = b.YSize
        ul = float(b.ReadAsArray(0, 0, 1, 1))
        #ur = float(b.ReadAsArray(ns-1, 0, 1, 1))
        lr = float(b.ReadAsArray(ns-1, nl-1, 1, 1))
        #ll = float(b.ReadAsArray(0, nl-1, 1, 1))
        #Probably better to use 3/4 corner criterion
        #if ul == ur == lr == ll:
        if np.isnan(ul) or ul == lr:
            b_ndv = ul
        else:
            #Assume ndv is 0
            b_ndv = 0
    elif np.isnan(b_ndv):
        b_dt = gdal.GetDataTypeName(b.DataType)
        if 'Float' in b_dt:
            b_ndv = np.nan
        else:
            b_ndv = 0
    return b_ndv

#Write out a recarray as a csv

def __init__(self,ds_filename,load_data=True,latlon=True,band=1,
        spatial_ref=None,geo_transform=None,downsampl=1):
        """ Construct object with raster from a single band dataset. 

        Parameters:
            ds_filename : filename of the dataset to import
            load_data : - True, to import the data into obj.r. 
                        - False, to not load any data.
                        - tuple (left, right, bottom, top) to load subset; 
                          obj.extent will be set to reflect subset area.
            latlon : default True. Only used if load_data=tuple. Set as False
                     if tuple is projected coordinates, True if WGS84.
            band : default 1. Specify GDAL band number to load. If you want to
                   load multiple bands at once use MultiBandRaster instead.
            downsampl : default 1. Used to down-sample the image when loading it. 
                A value of 2 for example will multiply the resolution by 2. 

            Optionally, you can manually specify/override the georeferencing. 
            To do this you must set both of the following parameters:

            spatial_ref : a OSR SpatialReference instance
            geo_transform : a Geographic Transform tuple of the form 
                            (top left x, w-e cell size, 0, top left y, 0, 
                             n-s cell size (-ve))

        """

        # Do basic dataset loading - set up georeferencing
        self._load_ds(ds_filename,spatial_ref=spatial_ref,
                      geo_transform=geo_transform)

        # Import band datatype
        band_tmp = self.ds.GetRasterBand(band)
        self.dtype = gdal.GetDataTypeName(band_tmp.DataType)

        # Load entire image
        if load_data == True:
            self.r = self.read_single_band(band,downsampl=downsampl)

        # Or load just a subset region
        elif isinstance(load_data,tuple) or isinstance(load_data,list):
            if len(load_data) == 4:
                self.r = self.read_single_band_subset(load_data,latlon=latlon,
                    band=band,update_info=True,downsampl=downsampl)

        elif load_data == False:
            return

        else:
            print('Warning : load_data argument not understood. No data loaded.')