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

def write_int16_to_tiff(name, data, sr, geot, nodata_val=None):
    assert data.dtype == np.int16
    gtiff_drv = gdal.GetDriverByName("GTiff")
    tiff_file = gtiff_drv.Create(name, data.shape[1], data.shape[0], 1,
                                 gdal.GDT_Int16,
                                 options=['COMPRESS=DEFLATE', 'ZLEVEL=1'])
    tiff_file.SetGeoTransform(geot)
    tiff_file.SetProjection(sr)

    band = tiff_file.GetRasterBand(1)
    if nodata_val is not None:
        band.SetNoDataValue(nodata_val)
    band.WriteArray(data)
    band.FlushCache()
    del band
    del tiff_file

def open_data(filename):
    '''
    The function open and load the image given its name. 
    The type of the data is checked from the file and the scipy array is initialized accordingly.
    Input:
        filename: the name of the file
    Output:
        im: the data cube
        GeoTransform: the geotransform information 
        Projection: the projection information
    '''
    data = gdal.Open(filename,gdal.GA_ReadOnly)
    if data is None:
        print 'Impossible to open '+filename
        exit()
    nc = data.RasterXSize
    nl = data.RasterYSize
    d  = data.RasterCount

    # Get the type of the data
    gdal_dt = data.GetRasterBand(1).DataType
    if gdal_dt == gdal.GDT_Byte:
        dt = 'uint8'
    elif gdal_dt == gdal.GDT_Int16:
        dt = 'int16'
    elif gdal_dt == gdal.GDT_UInt16:
        dt = 'uint16'
    elif gdal_dt == gdal.GDT_Int32:
        dt = 'int32'
    elif gdal_dt == gdal.GDT_UInt32:
        dt = 'uint32'

    elif gdal_dt == gdal.GDT_Float32:
        dt = 'float32'
    elif gdal_dt == gdal.GDT_Float64:
        dt = 'float64'
    elif gdal_dt == gdal.GDT_CInt16 or gdal_dt == gdal.GDT_CInt32 or gdal_dt == gdal.GDT_CFloat32 or gdal_dt == gdal.GDT_CFloat64 :
        dt = 'complex64'
    else:
        print 'Data type unkown'
        exit()

    # Initialize the array
    if d == 1:
        im = sp.empty((nl,nc),dtype=dt) 
    else:
        im = sp.empty((nl,nc,d),dtype=dt) 

    if d == 1:
        im[:,:]=data.GetRasterBand(1).ReadAsArray()
    else :
        for i in range(d):
            im[:,:,i]=data.GetRasterBand(i+1).ReadAsArray()

    GeoTransform = data.GetGeoTransform()
    Projection = data.GetProjection()
    data = None
    return im,GeoTransform,Projection

def write_data(outname,im,GeoTransform,Projection):
    '''
    The function write the image on the  hard drive.
    Input: 
        outname: the name of the file to be written
        im: the image cube
        GeoTransform: the geotransform information 
        Projection: the projection information
    Output:
        Nothing --
    '''
    nl = im.shape[0]
    nc = im.shape[1]
    if im.ndim == 2:
        d=1
    else:
        d = im.shape[2]

    driver = gdal.GetDriverByName('GTiff')
    dt = im.dtype.name
    # Get the data type
    if dt == 'bool' or dt == 'uint8':
        gdal_dt=gdal.GDT_Byte
    elif dt == 'int8' or dt == 'int16':
        gdal_dt=gdal.GDT_Int16
    elif dt == 'uint16':
        gdal_dt=gdal.GDT_UInt16
    elif dt == 'int32':
        gdal_dt=gdal.GDT_Int32
    elif dt == 'uint32':
        gdal_dt=gdal.GDT_UInt32
    elif dt == 'int64' or dt == 'uint64' or dt == 'float16' or dt == 'float32':
        gdal_dt=gdal.GDT_Float32
    elif dt == 'float64':
        gdal_dt=gdal.GDT_Float64
    elif dt == 'complex64':
        gdal_dt=gdal.GDT_CFloat64
    else:
        print 'Data type non-suported'
        exit()

    dst_ds = driver.Create(outname,nc,nl, d, gdal_dt)
    dst_ds.SetGeoTransform(GeoTransform)
    dst_ds.SetProjection(Projection)

    if d==1:
        out = dst_ds.GetRasterBand(1)
        out.WriteArray(im)
        out.FlushCache()
    else:
        for i in range(d):
            out = dst_ds.GetRasterBand(i+1)
            out.WriteArray(im[:,:,i])
            out.FlushCache()
    dst_ds = None

def warp(nimrod_dataset):
    """
    Warps the data array into one that has longitude/latitude as axes an fits
    the EPSG:4326 spatial reference system. The original array has the srs
    EPSG:27700 (OSGB 1936).
    :param nimrod_dataset: dictionary containing the data from the NIMROD file
    """

    # http://gis.stackexchange.com/questions/139906/replicating-result-of-gdalwarp-using-gdal-python-bindings
    # Create synthetic data
    gtiff_drv = gdal.GetDriverByName('MEM')
    cols, rows = nimrod_dataset['cols'], nimrod_dataset['rows']
    raster = np.reshape(nimrod_dataset['data'], (cols, rows))
    raster = np.int16(raster)
    top_left = (nimrod_dataset['start_easting'], nimrod_dataset['start_northing'])
    pixel_height = nimrod_dataset['column_interval']
    pixel_width = nimrod_dataset['row_interval']

    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(27700)

    src_geotran = [top_left[0], pixel_width,           0,
                   top_left[1],     0,      -pixel_height]

    rows, cols = raster.shape
    src_ds = gtiff_drv.Create(
        'test_epsg3413.tif',
        cols, rows, 1,
        gdal.GDT_Int16)
    src_ds.SetGeoTransform(src_geotran)
    src_ds.SetProjection(src_srs.ExportToWkt())
    src_ds.GetRasterBand(1).WriteArray(raster)

    # Transform to EPSG: 4326
    dest_srs = osr.SpatialReference()
    dest_srs.ImportFromEPSG(4326)

    int_ds = gdal.AutoCreateWarpedVRT(src_ds, src_srs.ExportToWkt(), dest_srs.ExportToWkt())
    nimrod_dataset['data_warped'] = int_ds.GetRasterBand(1).ReadAsArray()
    nimrod_dataset['GeoTransform'] = int_ds.GetGeoTransform()

    src_ds = None
    int_ds = None

    return nimrod_dataset