Python gdal.UseExceptions() Examples

def singleBand(self, rasterFp):
        gdal.UseExceptions()        
        '''打开栅格数据'''
        try:
            src_ds=gdal.Open(rasterFp)
        except RuntimeError as e:
            print( 'Unable to open %s'% rasterFp)
            sys.exit(1)
        #获取栅格信息
        rasterInfo={"RasterXSize":src_ds.RasterXSize,
                    "RasterYSize":src_ds.RasterYSize,
                    "RasterProjection":src_ds.GetProjection(),
                    "GeoTransform":src_ds.GetGeoTransform()}
        
        '''获取单波段像元值'''
        bandValue=src_ds.GetRasterBand(1).ReadAsArray().astype(np.float)
        print("readed rasterDate!")
        return bandValue,rasterInfo #返回该波段，为数组形式

#影像数据裁切 

def singleBand(self, rasterFp):
        gdal.UseExceptions()        
        '''打开栅格数据'''
        try:
            src_ds=gdal.Open(rasterFp)
        except RuntimeError as e:
            print( 'Unable to open %s'% rasterFp)
            sys.exit(1)
        #获取栅格信息
        rasterInfo={"RasterXSize":src_ds.RasterXSize,
                    "RasterYSize":src_ds.RasterYSize,
                    "RasterProjection":src_ds.GetProjection(),
                    "GeoTransform":src_ds.GetGeoTransform()}
        
        '''获取单波段像元值'''
        bandValue=src_ds.GetRasterBand(1).ReadAsArray().astype(np.float)
        print("readed rasterDate!")
        return bandValue,rasterInfo #返回该波段，为数组形式 

def singleBand(rasterFp):
    gdal.UseExceptions()        
    '''打开栅格数据'''
    try:
        src_ds=gdal.Open(rasterFp)
    except RuntimeError as e:
        # print( 'Unable to open %s'% rasterFp)
        sys.exit(1)
    #获取栅格信息
    rasterInfo={"RasterXSize":src_ds.RasterXSize,
                "RasterYSize":src_ds.RasterYSize,
                "RasterProjection":src_ds.GetProjection(),
                "GeoTransform":src_ds.GetGeoTransform()}
    
    '''获取单波段像元值'''
    bandValue=src_ds.GetRasterBand(1).ReadAsArray().astype(np.float)
    # print("readed rasterDate!")
    return bandValue,rasterInfo #返回该波段，为数组形式 

def rasterRW(self, LSTValue,resultsPath,LSTSavingFn,para):
        gdal.UseExceptions()    
    #    '''打开栅格数据'''
    #    try:
    #        src_ds=gdal.Open(os.path.join(resultsPath,LSTSavingFn))
    #    except RuntimeError as e:
    #        print( 'Unable to open %s'% os.path.join(resultsPath,LSTSavingFn))
    #        print(e)
    #        sys.exit(1)
    #    print("metadata:",src_ds.GetMetadata())   
      
        '''初始化输出栅格'''
        driver=gdal.GetDriverByName('GTiff')
        print(para['RasterXSize'],para['RasterYSize'])
        out_raster=driver.Create(os.path.join(resultsPath,LSTSavingFn),para['RasterXSize'],para['RasterYSize'],1,gdal.GDT_Float64)
        out_raster.SetProjection(para['RasterProjection']) #设置投影与参考栅格同
        out_raster.SetGeoTransform(para['GeoTransform']) #配置地理转换与参考栅格同
        
        '''将数组传给栅格波段，为栅格值'''
        out_band=out_raster.GetRasterBand(1)
        out_band.WriteArray(LSTValue)
        
    #    '''设置overview'''
    #    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
    #    out_raster.BuildOverviews('average', overviews)
        
        '''清理缓存与移除数据源'''
        out_band.FlushCache()
        out_band.ComputeStatistics(False)
    #    del src_ds,out_raster,out_band        
        del out_raster,out_band

##解译精度评价。采样的数据是使用GIS平台人工判断提取，样例文件在Github中获取。 

def singleRaster(fn,raster_lyr):
    gdal.UseExceptions()
    
    '''打开栅格数据'''
    try:
        src_ds=gdal.Open(os.path.join(fn,raster_lyr))
    except RuntimeError as e:
        print( 'Unable to open %s'% os.path.join(fn,raster_lyr))
        sys.exit(1)
    print("metadata:",src_ds.GetMetadata())       
    
    '''获取所有波段'''
    band=src_ds.GetRasterBand(1).ReadAsArray().astype(np.float)
    return band #返回该波段，为数组形式 

def loadModelPredict(studyRegionFn,modelPath,mask):
    gdal.UseExceptions()    
    '''打开栅格数据'''
    try:
        src_ds=gdal.Open(studyRegionFn)
    except RuntimeError as e:
        print( 'Unable to open %s'% studyRegionFn)
#        print(e)
        sys.exit(1)
    print("metadata:",src_ds.GetMetadata())       
    
    '''获取波段'''
    studyRegion=src_ds.GetRasterBand(1).ReadAsArray().astype(np.float)
 
    '''按照mask划分解释变量和目标变量'''
    image_shape=(row, col)    
    dataFlatten=np.copy(studyRegion).reshape(-1)  
    maskFlatten=mask.reshape(-1)
    maskBool=maskFlatten==1
    X=dataFlatten[maskBool].reshape(1,-1)
    
    model=joblib.load(modelPath) #加载已经训练好的回归模型，用于预测
    y_p=model.predict(X)
    
    print(X.std(),y_p.std(),X.max(),y_p.max())
    y_p_scaled=(y_p - y_p.min()) / (y_p.max() - y_p.min()) * (X.max() - X.min()) + X.min() #调整预测数据的区间与解释变量的区间同，也可忽略
    print(dataFlatten.shape,maskBool.shape,y_p_scaled.shape)
    
    true_pred=dataFlatten
    true_pred[maskBool==False]=y_p_scaled.reshape(-1)   #拼合解释变量和预测数据

    '''打印显示结果'''
    fig, (ax1, ax2) = plt.subplots(figsize=(17, 17), ncols=2)
    pre=ax1.imshow(studyRegion, cmap=plt.cm.RdYlGn, interpolation='none')
#    fig.colorbar(pre, ax=ax1)
    ori=ax2.imshow(true_pred.reshape(image_shape), cmap=plt.cm.RdYlGn, interpolation='none')
#    fig.colorbar(ori, ax=ax2)
    plt.show() 
    
    return true_pred.reshape(image_shape) 

def rasterRW(self, rasterValue,resultsPath,rasterSavingFn,para):
        gdal.UseExceptions()    
    #    '''打开栅格数据'''
    #    try:
    #        src_ds=gdal.Open(os.path.join(resultsPath,rasterSavingFn))
    #    except RuntimeError as e:
    #        print( 'Unable to open %s'% os.path.join(resultsPath,rasterSavingFn))
    #        print(e)
    #        sys.exit(1)
    #    print("metadata:",src_ds.GetMetadata())   
      
        '''初始化输出栅格'''
        driver=gdal.GetDriverByName('GTiff')
        print(para['RasterXSize'],para['RasterYSize'])
        out_raster=driver.Create(os.path.join(resultsPath,rasterSavingFn),para['RasterXSize'],para['RasterYSize'],1,gdal.GDT_Float64)
        out_raster.SetProjection(para['RasterProjection']) #设置投影与参考栅格同
        out_raster.SetGeoTransform(para['GeoTransform']) #配置地理转换与参考栅格同
        
        '''将数组传给栅格波段，为栅格值'''
        out_band=out_raster.GetRasterBand(1)
        out_band.WriteArray(rasterValue)
        
    #    '''设置overview'''
    #    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
    #    out_raster.BuildOverviews('average', overviews)
        
        '''清理缓存与移除数据源'''
        out_band.FlushCache()
        out_band.ComputeStatistics(False)
    #    del src_ds,out_raster,out_band        
        del out_raster,out_band 

def rasterRW(rasterValue,resultsPath,rasterSavingFn,para):
    gdal.UseExceptions()    
#    '''打开栅格数据'''
#    try:
#        src_ds=gdal.Open(os.path.join(resultsPath,rasterSavingFn))
#    except RuntimeError as e:
#        print( 'Unable to open %s'% os.path.join(resultsPath,rasterSavingFn))
#        print(e)
#        sys.exit(1)
#    print("metadata:",src_ds.GetMetadata())   
  
    '''初始化输出栅格'''
    driver=gdal.GetDriverByName('GTiff')
    # print(para['RasterXSize'],para['RasterYSize'])
    out_raster=driver.Create(os.path.join(resultsPath,rasterSavingFn),para['RasterXSize'],para['RasterYSize'],1,gdal.GDT_Float64)
    out_raster.SetProjection(para['RasterProjection']) #设置投影与参考栅格同
    out_raster.SetGeoTransform(para['GeoTransform']) #配置地理转换与参考栅格同
    
    '''将数组传给栅格波段，为栅格值'''
    out_band=out_raster.GetRasterBand(1)
    out_band.WriteArray(rasterValue)
    
#    '''设置overview'''
#    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
#    out_raster.BuildOverviews('average', overviews)
    
    '''清理缓存与移除数据源'''
    out_band.FlushCache()
    out_band.ComputeStatistics(False)
#    del src_ds,out_raster,out_band        
    del out_raster,out_band 

def rasterRW(self, rasterArray,resultsPath,resultsFn,para):
        gdal.UseExceptions()    
    #    '''打开栅格数据'''
    #    try:
    #        src_ds=gdal.Open(os.path.join(resultsPath,resultsFn))
    #    except RuntimeError as e:
    #        print( 'Unable to open %s'% os.path.join(resultsPath,resultsFn))
    #        print(e)
    #        sys.exit(1)
    #    print("metadata:",src_ds.GetMetadata())   
      
        '''初始化输出栅格'''
        driver=gdal.GetDriverByName('GTiff')
        print(para['RasterXSize'],para['RasterYSize'])
        out_raster=driver.Create(os.path.join(resultsPath,resultsFn),para['RasterXSize'],para['RasterYSize'],1,gdal.GDT_Float64)
        out_raster.SetProjection(para['RasterProjection']) #设置投影与参考栅格同
        out_raster.SetGeoTransform(para['GeoTransform']) #配置地理转换与参考栅格同
        
        '''将数组传给栅格波段，为栅格值'''
        out_band=out_raster.GetRasterBand(1)
        out_band.WriteArray(rasterArray)
        
    #    '''设置overview'''
    #    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
    #    out_raster.BuildOverviews('average', overviews)
        
        '''清理缓存与移除数据源'''
        out_band.FlushCache()
        out_band.ComputeStatistics(False)
    #    del src_ds,out_raster,out_band    
        print("raster saved successfully!")    
        del out_raster,out_band 

def export_geomap2geotiff(path, geo_model, geo_map=None, geotiff_filepath=None):
    """

    Args:
        path (str): Filepath for the exported geotiff, must end in .tif
        geo_map (np.ndarray): 2-D array containing the geological map
        cmap (matplotlib colormap): The colormap to be used for the export
        geotiff_filepath (str): Filepath of the template geotiff

    Returns:
        Saves the geological map as a geotiff to the given path.
    """
    import gdal

    plot = PlotData2D(geo_model)
    cmap = plot._cmap
    norm = plot._norm

    if geo_map is None:
        geo_map = geo_model.solutions.geological_map[0].reshape(geo_model._grid.topography.resolution)

    if geotiff_filepath is None:
        # call the other function
        print('stupid')

    # **********************************************************************
    geo_map_rgb = SM(norm=norm, cmap=cmap).to_rgba(geo_map.T) # r,g,b,alpha
    # **********************************************************************
    # gdal.UseExceptions()
    ds = gdal.Open(geotiff_filepath)
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray()
    [cols, rows] = arr.shape

    outFileName = path
    driver = gdal.GetDriverByName("GTiff")
    options = ['PROFILE=GeoTiff', 'PHOTOMETRIC=RGB', 'COMPRESS=JPEG']
    outdata = driver.Create(outFileName, rows, cols, 3, gdal.GDT_Byte, options=options)

    outdata.SetGeoTransform(ds.GetGeoTransform())  # sets same geotransform as input
    outdata.SetProjection(ds.GetProjection())  # sets same projection as input
    outdata.GetRasterBand(1).WriteArray(geo_map_rgb[:, ::-1, 0].T * 256)
    outdata.GetRasterBand(2).WriteArray(geo_map_rgb[:, ::-1, 1].T * 256)
    outdata.GetRasterBand(3).WriteArray(geo_map_rgb[:, ::-1, 2].T * 256)
    outdata.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
    outdata.GetRasterBand(2).SetColorInterpretation(gdal.GCI_GreenBand)
    outdata.GetRasterBand(3).SetColorInterpretation(gdal.GCI_BlueBand)
    # outdata.GetRasterBand(4).SetColorInterpretation(gdal.GCI_AlphaBand)  # alpha band

    # outdata.GetRasterBand(1).SetNoDataValue(999)##if you want these values transparent
    outdata.FlushCache()  # saves to disk
    outdata = None  # closes file (important)
    band = None
    ds = None

    print("Successfully exported geological map to  " +path) 

def rasterCal(fn,raster_lyr,raster_lyr_w):   
    gdal.UseExceptions()
    
    '''打开栅格数据'''
    try:
        src_ds=gdal.Open(os.path.join(fn,raster_lyr))
    except RuntimeError as e:
        print( 'Unable to open %s'% os.path.join(fn,raster_lyr))
        print(e)
        sys.exit(1)
    print("metadata:",src_ds.GetMetadata())   
    
    '''获取所有波段'''
    srcband=[]
    for band_num in range(1,8):
        try:
            srcband.append(src_ds.GetRasterBand(band_num))
        except RuntimeError as e:
            print('Band ( %i ) not found' % band_num)
            print(e)
            sys.exit(1)
    print(srcband)
    
    '''获取用于NDVI计算的红和近红波段数组,并计算ndvi'''
    red=srcband[3].ReadAsArray().astype(np.float)
    nir=srcband[4].ReadAsArray()
    red=np.ma.masked_where(nir+red==0,red)  #确定分母不为零
    ndvi=(nir-red)/(nir+red)
    ndvi=ndvi.filled(-99)
    print(ndvi.shape,ndvi.std(),ndvi.max(),ndvi.min(),ndvi.mean())
    
    '''初始化输出栅格'''
    driver=gdal.GetDriverByName('GTiff')
    if os.path.exists(os.path.join(fn,raster_lyr_w)):
        driver.Delete(os.path.join(fn,raster_lyr_w))
    out_raster=driver.Create(os.path.join(fn,raster_lyr_w),src_ds.RasterXSize,src_ds.RasterYSize,1,gdal.GDT_Float64)
    out_raster.SetProjection(src_ds.GetProjection()) #设置投影与参考栅格同
    out_raster.SetGeoTransform(src_ds.GetGeoTransform()) #配置地理转换与参考栅格同
    
    '''将数组传给栅格波段，为栅格值'''
    out_band=out_raster.GetRasterBand(1)
    out_band.WriteArray(ndvi)
    
    '''设置overview'''
    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
    out_raster.BuildOverviews('average', overviews)
    
    '''清理缓存与移除数据源'''
    out_band.FlushCache()
    out_band.ComputeStatistics(False)
    del src_ds,out_raster,out_band
    return ndvi 

def rasterRW(fn,raster_lyr,raster_lyr_w):
    gdal.UseExceptions()
    
    '''打开栅格数据'''
    try:
        src_ds=gdal.Open(os.path.join(fn,raster_lyr))
    except RuntimeError as e:
        print( 'Unable to open %s'% os.path.join(fn,raster_lyr))
        print(e)
        sys.exit(1)
    print("metadata:",src_ds.GetMetadata())   
    
    '''获取所有波段'''
    srcband=[]
    for band_num in range(1,5):
        try:
            srcband.append(src_ds.GetRasterBand(band_num))
        except RuntimeError as e:
            print('Band ( %i ) not found' % band_num)
            print(e)
            sys.exit(1)
    print(srcband)
    
    '''获取用于NDVI计算的红和近红波段数组,并计算ndvi'''
    red=srcband[0].ReadAsArray().astype(np.float)
    nir=srcband[3].ReadAsArray()
    red=np.ma.masked_where(nir+red==0,red)
    ndvi=(nir-red)/(nir+red)
    ndvi=ndvi.filled(-99)
    print(ndvi.shape,ndvi.std())
    
    '''初始化输出栅格'''
    driver=gdal.GetDriverByName('GTiff')
    out_raster=driver.Create(os.path.join(fn,raster_lyr_w),src_ds.RasterXSize,src_ds.RasterYSize,1,gdal.GDT_Float64)
    out_raster.SetProjection(src_ds.GetProjection()) #设置投影与参考栅格同
    out_raster.SetGeoTransform(src_ds.GetGeoTransform()) #配置地理转换与参考栅格同
    
    '''将数组传给栅格波段，为栅格值'''
    out_band=out_raster.GetRasterBand(1)
    out_band.WriteArray(ndvi)
    
    '''设置overview'''
    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
    out_raster.BuildOverviews('average', overviews)
    
    '''清理缓存与移除数据源'''
    out_band.FlushCache()
    out_band.ComputeStatistics(False)
    del src_ds,out_raster,out_band