Python ogr.Open() Examples

def rasterize(data, vectorSrc, field, outFile):
    dataSrc = gdal.Open(data)
    import ogr
    shp = ogr.Open(vectorSrc)

    lyr = shp.GetLayer()

    driver = gdal.GetDriverByName('GTiff')
    dst_ds = driver.Create(
        outFile,
        dataSrc.RasterXSize,
        dataSrc.RasterYSize,
        1,
        gdal.GDT_UInt16)
    dst_ds.SetGeoTransform(dataSrc.GetGeoTransform())
    dst_ds.SetProjection(dataSrc.GetProjection())
    if field is None:
        gdal.RasterizeLayer(dst_ds, [1], lyr, None)
    else:
        OPTIONS = ['ATTRIBUTE=' + field]
        gdal.RasterizeLayer(dst_ds, [1], lyr, None, options=OPTIONS)

    data, dst_ds, shp, lyr = None, None, None, None
    return outFile 

def pointReading(fn,pt_lyrName_r):
    ds=ogr.Open(fn,0) #0为只读模式，1为编辑模式
    if ds is None:
        sys.exit('Could not open{0}'.format(fn))
    pt_lyr=ds.GetLayer(pt_lyrName_r) #可以直接数据层(文件)名或者指定索引
    vp = VectorPlotter(True)  #显示vector数据
    vp.plot(pt_lyr,'bo')
    
    i=0
    for feat in pt_lyr: #循环feature
        pt=feat.geometry()
        pt_x=pt.GetX()
        pt_y=pt.GetY()
        name=feat.GetField('NAME')
        kind=feat.GetField('KIND')
        print(name,kind,(pt_x,pt_y,))
        i+=1
        if i==12:
            break
    del ds 

def reproject_modis_to_geotiff(input_folder, dest_prj=32643):
    '''
    Function to convert all modis hdf in folder to geotiff using gdal
    required libraries: osr, gdal
    :param input_folder: where hds files are stored
    :param dest_prj: default is wgs 84 / UTM 43N (EPSG 32643)
    '''
    files_list = os.listdir(input_folder)
    for item in files_list:
        if fnmatch.fnmatch(item, '*.hdf'):
            input_file_name = input_folder + '/' + item
            output_file_name = input_folder + '/' + item[0:23]
            img = gdal.Open(input_file_name)
            subset = img.GetSubDatasets()
            in_raster = subset[0][0]
            reproject_dataset(dataset=in_raster, output_file_name=output_file_name+ '.tif')

# reproject
# reproject_modis_to_geotiff(input_folder="/media/kiruba/New Volume/MODIS/ET/scratch")

# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL. 

def test_stack_images(managed_multiple_geotiff_dir):
    test_dir = managed_multiple_geotiff_dir
    images = [os.path.join(test_dir.path, image) for image in os.listdir(test_dir.path)]
    images.sort()
    result_path = os.path.join(test_dir.path, "test_out.tif") # This isn't turning up in the expected order
    pyeo.raster_manipulation.stack_images(images, result_path)
    result = gdal.Open(result_path)
    assert result.ReadAsArray().shape[0] == 4  # If result has 4 bands, success
    assert result.ReadAsArray()[0,0,4] == 3
    assert result.ReadAsArray()[1,0,4] == 13 

def polygonReading(fn,pg_lyrName_r):    
    ds=ogr.Open(fn,0) #0为只读模式，1为编辑模式
    if ds is None:
        sys.exit('Could not open{0}'.format(fn))
    pg_lyr=ds.GetLayer(pg_lyrName_r) #可以直接数据层(文件)名或者指定索引
    
    vp = VectorPlotter(True) #显示vector数据
    vp.plot(pg_lyr)
    
    pg_schema=pg_lyr.schema #查看属性表字段名和类型
    for field in pg_schema:
        print(field.name,field.GetTypeName())     
        
    i=0    
    for feat in pg_lyr: #循环feature
        print("..................................................................")
        atts=feat.items()
        print(atts)
        pg=feat.geometry() #获取feature的几何对象
#        ring=
        name=feat.GetField('NAME')  #读取feature的属性
        Shape_Area=feat.GetField('Shape_Area')
        print(name,Shape_Area,'\n',pg,'\n')   
        
        for j in range(pg.GetGeometryCount()):  #循环几何对象获取ring,获取顶点坐标
            ring=pg.GetGeometryRef(j)
            for coordi in ring.GetPoints():
                print(coordi)
        i+=1
        if i==12:
            break
    del ds 

def lineReading(fn,ln_lyrName_r):    
    ds=ogr.Open(fn,0) #0为只读模式，1为编辑模式
    if ds is None:
        sys.exit('Could not open{0}'.format(fn))
    ln_lyr=ds.GetLayer(ln_lyrName_r) #可以直接数据层(文件)名或者指定索引
    
    vp = VectorPlotter(True) #显示vector数据
    vp.plot(ln_lyr,'bo')
    
    ln_schema=ln_lyr.schema #查看属性表字段名和类型
    for field in ln_schema:
        print(field.name,field.GetTypeName())     
        
    i=0    
    for feat in ln_lyr: #循环feature
        ln=feat.geometry() #获取feature的几何对象
        name=feat.GetField('name')  #读取feature的属性
        Shape_Leng=feat.GetField('Shape_Leng')
        print(name,Shape_Leng,'\n',ln,'\n',ln.GetPointCount())        
        for j in range(ln.GetPointCount()):  #循环几何对象(线)d的vertex顶点
            if j<6:
                print((i,ln.GetX(i),ln.GetY(i))) #只能通过GetX()和GetY()的方法获取顶点坐标
        i+=1
        if i==12:
            break
    del ds 

def gettingMetadata(fn,lyrName):
    ds=ogr.Open(fn,0)
    if ds is None:
        sys.exit('Could not open{0}'.format(fn))
    lyr=ds.GetLayer(lyrName)    
    print("extent:",lyr.GetExtent()) #(min_x,max_x, min_y, max_y)
    print("GeomType:",lyr.GetGeomType(),"wkbPoint?",lyr.GetGeomType()==ogr.wkbPoint) #返回索引，根据索引查找对应的数据类型
    feat=lyr.GetFeature(0)
    print("GeomTypeByfeature:",feat.geometry().GetGeometryName())
    print("spatialRef:",lyr.GetSpatialRef())  #空间坐标
    print("fieldAttri:",[(field.name,field.GetTypeName()) for field in lyr.schema])
    lyrName=lyrName+".shp"
    pb.print_attributes(os.path.join(fn,lyrName), 3, ["NAME","KIND"]) 

def pts2raster(shapefile,RASTER_PATH,cellSize,field_name=False):
    from osgeo import gdal, ogr

    # Define pixel_size and NoData value of new raster
    pixel_size = cellSize
    NoData_value = -9999
    
    # Filename of input OGR file
    vector_ptsShp_fn = shapefile
    
    # Filename of the raster Tiff that will be created
    raster_ptsShp_fn = RASTER_PATH
    
    # Open the data source and read in the extent
    source_ds = ogr.Open(vector_ptsShp_fn)
    source_layer = source_ds.GetLayer()
    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    
    # Create the destination data source
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_ptsShp_fn, x_res, y_res, 1, gdal.GDT_Int32 )
    target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(NoData_value)
    
    # Rasterize
    # gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])
    # Rasterize
    if field_name:
        gdal.RasterizeLayer(target_ds,[1], source_layer,options=["ATTRIBUTE={0}".format(field_name)])
        # print("write:",field_name)
    else:
        gdal.RasterizeLayer(target_ds,[1], source_layer,burn_values=[-1])   
    return gdal.Open(RASTER_PATH).ReadAsArray() 

#批量计算 

def _get_gdb_ds(self):
        """Get the geodatabase OGR data source object."""
        if not self.arcpy_found:
            return ogr.Open(self.path, 0)

    # ---------------------------------------------------------------------- 

def __init__(self, source_file):
        self._ds = ogr.Open(source_file)
        self._targetGeomColSRID = -1
        self._geomType = ''
        self._dbSession = STDMDb.instance().session
        self._mapped_cls = None
        self._mapped_doc_cls = None
        self._current_profile = current_profile()
        self._source_doc_manager = None 

def shpGeom2pg(self, aoi, tablename):
        """
        This function is a wrapper to import a shapefile geometry to a PostGreSQL database
        """

        sqlCmd = 'DROP TABLE IF EXISTS {}'.format(tablename)
        self.cursor.execute(sqlCmd)

        fList = 'id smallint, geometry geometry'
        sqlCmd = 'CREATE TABLE {} ({})'.format(tablename, fList)
        self.cursor.execute(sqlCmd)

        prjFile = '{}.prj'.format(aoi[:-4])
        inProj4 = get_proj4(prjFile)

        sf = ogr.Open(aoi)
        layer = sf.GetLayer(0)
        for i in range(layer.GetFeatureCount()):
            feature = layer.GetFeature(i)
            wkt = feature.GetGeometryRef().ExportToWkt()

            if inProj4 is not '+proj=longlat +datum=WGS84 +no_defs':
                wkt = reproject_geometry(wkt, inProj4, 4326)

            wkt = 'St_GeomFromText(\'{}\', 4326)'.format(wkt)
            values = '(\'{}\', {})'.format(i, wkt)
            sql_cmd = 'INSERT INTO {} VALUES {}'.format(tablename, values)
            self.cursor.execute(sql_cmd) 

def kml_to_wkt(kmlfile):

    shp = ogr.Open(os.path.abspath(kmlfile))
    lyr = shp.GetLayerByName()
    for feat in lyr:
        geom = feat.GetGeometryRef()
    wkt = str(geom)

    return wkt 

def shp_to_wkt(shapefile, buffer=None, convex=False, envelope=False):
    '''A helper function to translate a shapefile into WKT


    '''

    # get filepaths and proj4 string
    shpfile = os.path.abspath(shapefile)
    prjfile = shpfile[:-4] + '.prj'
    proj4 = get_proj4(prjfile)

    lyr_name = os.path.basename(shapefile)[:-4]
    shp = ogr.Open(os.path.abspath(shapefile))
    lyr = shp.GetLayerByName(lyr_name)
    geom = ogr.Geometry(ogr.wkbGeometryCollection)

    for feat in lyr:
        geom.AddGeometry(feat.GetGeometryRef())
        wkt = geom.ExportToWkt()

    if proj4 != '+proj=longlat +datum=WGS84 +no_defs':
        print(' INFO: Reprojecting AOI file to Lat/Long (WGS84)')
        wkt = reproject_geometry(wkt, proj4, 4326).ExportToWkt()

    # do manipulations if needed
    wkt = wkt_manipulations(wkt, buffer=buffer, convex=convex,
                            envelope=envelope)

    return wkt 

def get_proj4(prjfile):
    '''Get the proj4 string from a projection file of a shapefile

    Args:
        prjfile: a .prj file of a shapefile

    Returns:
        str: PROJ4 code

    '''

    prj_file = open(prjfile, 'r')
    prj_string = prj_file.read()

    # Lambert error
    if '\"Lambert_Conformal_Conic\"' in prj_string:

        print(' ERROR: It seems you used an ESRI generated shapefile'
              ' with Lambert Conformal Conic projection. ')
        print(' This one is not compatible with Open Standard OGR/GDAL'
              ' tools used here. ')
        print(' Reproject your shapefile to a standard Lat/Long projection'
              ' and try again')
        exit(1)

    srs = osr.SpatialReference()
    srs.ImportFromESRI([prj_string])
    return srs.ExportToProj4() 

def __init__(self, inShape, inField, outValidation,
                 outTrain, number=50, percent=True):
        """
        inShape : str path file (e.g. '/doc/ref.shp')
        inField : string column name (e.g. 'class')
        outValidation : str path of shp output file (e.g. '/tmp/valid.shp')
        outTrain : str path of shp output file (e.g. '/tmp/train.shp')
        """
        from sklearn.model_selection import train_test_split

        if percent:
            number = number / 100.0
        else:
            number = int(number)

        lyr = ogr.Open(inShape)
        lyr1 = lyr.GetLayer()
        FIDs = np.zeros(lyr1.GetFeatureCount(), dtype=int)
        Features = []
        #unselFeat = []
        #current = 0

        for i, j in enumerate(lyr1):
            #print j.GetField(inField)
            FIDs[i] = j.GetField(inField)
            Features.append(j)
            #current += 1
        srs = lyr1.GetSpatialRef()
        lyr1.ResetReading()

        ##
        if percent:
            validation, train = train_test_split(
                Features, test_size=number, train_size=1 - number, stratify=FIDs)
        else:
            validation, train = train_test_split(
                Features, test_size=number, stratify=FIDs)

        self.saveToShape(validation, srs, outValidation)
        self.saveToShape(train, srs, outTrain) 

def open_data_band(filename):
    """!@brief The function open and load the image given its name.
    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:
            data : the opened data with gdal.Open() method
            im : empty table with right dimension (array)

    """
    data = gdal.Open(filename, gdal.GA_Update)
    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
    dt = getDTfromGDAL(gdal_dt)

    # Initialize the array
    im = np.empty((nl, nc), dtype=dt)
    return data, im 

def process_modis(input_folder, output_folder, scale_factor=0.1, null_value=32760, file_extension='*.tif'):
    files_list = os.listdir(input_folder)
    for item in files_list:
        if fnmatch.fnmatch(item, file_extension):
            in_raster = input_folder + '/' + item
            in_raster = gdal.Open(in_raster, GA_ReadOnly)
            out_raster = output_folder + '/' + 'proc_' + item
            band1 = in_raster.GetRasterBand(1)
            in_array = BandReadAsArray(band1)
            in_array[in_array > null_value] = np.nan
            data_out = in_array * scale_factor
            gdalnumeric.SaveArray(data_out, filename=out_raster, format='GTiff', prototype=in_raster) 

def test_stack_and_trim_images(managed_noncontiguous_geotiff_dir):
    # Test data is two five band 11x12 pixel geotiffs and a 10x10 polygon
    # The geotiffs upper left corners ar at 90,90 and 100,100
    test_dir = managed_noncontiguous_geotiff_dir
    old_image_path = os.path.join(test_dir.path, "se_test")
    new_image_path = os.path.join(test_dir.path, "ne_test")
    aoi_path = os.path.join(test_dir.path, "aoi")
    result_path = os.path.join(test_dir.path, "test_out.tif")
    pyeo.raster_manipulation.stack_and_trim_images(old_image_path, new_image_path, aoi_path, result_path)
    result = gdal.Open(result_path).ReadAsArray()
    assert result.shape == (10, 10, 10)  # Result should have 10 bands and be 10 by 10 

def test_pixel_bounds_from_polygon(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    layer = aoi.GetLayer(0)
    aoi_feature = layer.GetFeature(0)
    polygon = aoi_feature.GetGeometryRef()
    result = pyeo.coordinate_manipulation.pixel_bounds_from_polygon(raster, polygon)
    assert result == (1, 11, 1, 11) 

def test_point_to_pixel_coordinates(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    point = r"POINT (101 112)"
    raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    out = pyeo.coordinate_manipulation.point_to_pixel_coordinates(raster, point)
    # Since ne_test is an 11 by 12 image, tl corner 0,0 w 10m pixels,
    # we'd expect the br coords to be 10, 11
    assert out == (10, 11) 

def test_check_overlap(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    assert pyeo.coordinate_manipulation.check_overlap(test_raster, test_aoi) is True 

def test_get_raster_bounds(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    target = (0, 110, 0, 120)
    result = pyeo.coordinate_manipulation.get_raster_bounds(test_raster)
    assert result.GetEnvelope() == target 

def test_get_aoi_bounds(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    target = (10, 20, 10, 20)
    result = pyeo.coordinate_manipulation.get_aoi_bounds(test_aoi)
    assert result.GetEnvelope() == target 

def test_get_aoi_intersection(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "se_test"))
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    result = pyeo.coordinate_manipulation.get_aoi_intersection(test_raster, test_aoi)
    assert result.GetEnvelope() == (10, 20, 10, 20) 

def test_aoi_to_mask(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    result_path = os.path.join(test_dir.path, "test_out.tif")
    pyeo.aoi_to_mask(test_aoi, result_path)
    result = gdal.Open(result_path)
    assert result.ReadAsArray()[2, 2] == 1 

def clip_raster_by_vector(input_folder, mask_shapefile, output_folder, file_extension='*.tif', t_srs='EPSG:32643', no_data=32767 ):
    files_list = os.listdir(input_folder)
    ds = ogr.Open(mask_shapefile)
    lyr = ds.GetLayer(0)
    lyr.ResetReading()
    ft = lyr.GetNextFeature()
    for item in files_list:
        if fnmatch.fnmatch(item, file_extension):
            in_raster = input_folder + '/' + item
            out_raster = output_folder + '/' +'tg_' + item
            subprocess.call(['gdalwarp', in_raster, out_raster, '-cutline', mask_shapefile, '-t_srs', t_srs, '-crop_to_cutline', '-dstnodata', "%s" %no_data]) 

def reproject_dataset(dataset, output_file_name, wkt_from="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs", epsg_to=32643, pixel_spacing=1000, file_format='GTiff'):
    '''
    :param dataset: Modis subset name use gdal.GetSubdatasets()
    :param output_file_name: file location along with proper extension
    :param wkt_from: Modis wkt (default)
    :param epsg_to: integer default(32643)
    :param pixel_spacing: in metres
    :param file_format: default GTiff
    :return:
    '''
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_to)
    modis = osr.SpatialReference()
    modis.ImportFromProj4(wkt_from)
    tx = osr.CoordinateTransformation(modis, wgs84)
    g = gdal.Open(dataset)
    geo_t = g.GetGeoTransform()
    print geo_t
    x_size = g.RasterXSize
    y_size = g.RasterYSize
    (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
    (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + (geo_t[1]*x_size), geo_t[3]+ (geo_t[5]*y_size))
    mem_drv = gdal.GetDriverByName(file_format)
    dest = mem_drv.Create(output_file_name, int((lrx-ulx)/pixel_spacing), int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    new_geo = ([ulx, pixel_spacing, geo_t[2], uly, geo_t[4], -pixel_spacing])
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(wgs84.ExportToWkt())
    gdal.ReprojectImage(g, dest, modis.ExportToWkt(), wgs84.ExportToWkt(), gdal.GRA_Bilinear)
    print "reprojected" 

def clip_raster_by_vector(input_folder, mask_shapefile, output_folder, file_extension='*.tif', t_srs='EPSG:32643', no_data=32767 ):
    files_list = os.listdir(input_folder)
    ds = ogr.Open(mask_shapefile)
    lyr = ds.GetLayer(0)
    lyr.ResetReading()
    ft = lyr.GetNextFeature()
    for item in files_list:
        if fnmatch.fnmatch(item, file_extension):
            in_raster = input_folder + '/' + item
            out_raster = output_folder + '/' +'tg_' + item
            subprocess.call(['gdalwarp', in_raster, out_raster, '-cutline', mask_shapefile, '-t_srs', t_srs, '-crop_to_cutline', '-dstnodata', "%s" %no_data]) 

def get_feature_classes(self):
        """Get geodatabase feature classes as ordered dicts."""
        fcs = []
        if self.arcpy_found:
            arcpy.env.workspace = self.path
            # iterate feature classes within feature datasets
            fds = [fd for fd in arcpy.ListDatasets(feature_type='feature')]
            if fds:
                for fd in fds:
                    arcpy.env.workspace = os.path.join(self.path, fd)
                    for fc in arcpy.ListFeatureClasses():
                        od = self._get_fc_props(fc)
                        od['Feature dataset'] = fd
                        fcs.append(od)

            # iterate feature classes in the geodatabase root
            arcpy.env.workspace = self.path
            for fc in arcpy.ListFeatureClasses():
                od = self._get_fc_props(fc)
                fcs.append(od)

        else:
            ds = ogr.Open(self.path, 0)
            fcs_names = [
                ds.GetLayerByIndex(i).GetName()
                for i in range(0, ds.GetLayerCount())
                if ds.GetLayerByIndex(i).GetGeometryColumn()
            ]
            for fc_name in fcs_names:
                try:
                    fc_instance = FeatureClassOgr(self, fc_name)
                    od = OrderedDict()
                    for k, v in GDB_FC_PROPS.items():
                        od[v] = getattr(fc_instance, k, '')
                    # custom props
                    od['Row count'] = fc_instance.get_row_count()
                    fcs.append(od)
                except Exception as e:
                    print(e)
        return fcs

    # ---------------------------------------------------------------------- 

def dbf_to_df(path):
    """
    Return DataFrame from attributes stored in dBase/xBase format.

    Uses OGR's ESRI Shapefile driver to read records from the file.

    Parameters
    ----------
    path : str
        File path to the dBase/xBase file.

    Returns
    -------
    df : pandas.DataFrame

    """
    import ogr

    # Open the file and collect information on fields.
    dbf = ogr.Open(path)
    table = dbf.GetLayer()
    header = table.GetLayerDefn()
    ncolumns = header.GetFieldCount()
    column_names = [header.GetFieldDefn(i).GetName() for i in range(ncolumns)]
    column_types = [header.GetFieldDefn(i).GetType() for i in range(ncolumns)]

    def read(row, i):
        """Return i-th field of a record."""
        # For performance, use the appropriate field type function.
        fld_type = column_types[i]
        if fld_type == ogr.OFTInteger:
            return row.GetFieldAsInteger(i)
        elif fld_type == ogr.OFTReal:
            return row.GetFieldAsDouble(i)
        elif fld_type == ogr.OFTStringList:
            return row.GetFieldAsStringList(i)
        elif fld_type == ogr.OFTIntegerList:
            return row.GetFieldAsIntegerList(i)
        elif fld_type == ogr.OFTRealList:
            return row.GetFieldAsDoubleList(i)
        else:
            return row.GetFieldAsString(i)

    # Represent records with memory-efficient generators.
    values = lambda row: (read(row, i) for i in range(ncolumns))
    records = (values(row) for row in table)

    df = pd.DataFrame.from_records(records, columns=column_names,
                                   coerce_float=False)
    return df