常用的Python代码片段（地理相关）

把pandas的dataframe转为geopandas的地理格式

def df2gdf(df, lon_col='longitude', lat_col='latitude', epsg=4326, crs=None):
    gdf = gpd.GeoDataFrame(df)
    gdf['geometry'] = gpd.points_from_xy(x=df[lon_col], y=df[lat_col])
    gdf.geometry = gdf['geometry']
    if crs is None:
        gdf.set_crs(epsg=epsg, inplace=True)
    else:
        gdf.set_crs(crs=crs, inplace=True)
    return gdf

geopandas to pandas

df1 = pd.DataFrame(gdf.drop(columns='geometry'))
# for older versions of pandas (< 0.21), the drop part: gdf.drop('geometry', axis=1)

重采样

from osgeo import gdal
import glob
import os
os.chdir("/mnt/f/analysis/")
def resample(image, new_xres, new_yres, save_path):
    """
    Decrease the pixel size of the raster.
    Args:
        new_xres (int): desired resolution in x-direction
        new_yres (int): desired resolution in y-direction
        save_path (str): filepath to where the output file should be stored
    Returns: Nothing, it writes a raster file with decreased resolution.
    """
    props = image.GetGeoTransform()
    print('current pixel xsize:', props[1], 'current pixel ysize:', -props[-1])
    options = gdal.WarpOptions(options=['tr'], xRes=new_xres, yRes=new_yres)
    newfile = gdal.Warp(save_path, image, options=options)
    newprops = newfile.GetGeoTransform()
    print('new pixel xsize:', newprops[1], 'new pixel ysize:', -newprops[-1])
    print('file saved in ' + save_path)
IN_DIR = "CDL-clip"
OUT_DIR = "CDL_resample"
if not os.path.exists(OUT_DIR):
    os.mkdir(OUT_DIR)
for i in glob.glob(f"{IN_DIR}/*.tif"):
    print(f"[{i}]","* "*5)
    resample(gdal.Open(i), 0.008983152841195214,0.008983152841195214, i.replace(IN_DIR,OUT_DIR))

栅格运算

os.system(f"gdal_calc.py --calc \(B=={CROP_CLASS}\)*A --outfile {out_filename} \
        -A {IN_FILENAMES[i]} -B {MASK_FILENAMES[i]} \
        --extent union  --A_band 1 --B_band 1  --quiet --NoDataValue=0")

但是在windows下的gdal运行出错，在ubuntu系统内运行成功

读取地理栅格图像

def read_tif(tif_path, display=False):
    dataset = rio.open(tif_path)
    if display:
        print(f"width:{dataset.width}, height:{dataset.height}, bands:{dataset.count}")
      
        types = {i: dtype for i, dtype in zip(dataset.indexes, dataset.dtypes)}
        for k in types.keys():
            print(f"band_{k}: {types[k]}")
       
        print(f"{dataset.crs}")
        print(dataset.transform)
   
    return dataset
# 根据经纬度，从栅格中提取值（适用于EPSG:4326,其他坐标系未做测试）
def extract_value_from_tif(rs, lon,lat):
    '''
    geometry: geodataframe的geometry列
    rs: 栅格数据
    '''
    x = lon
    y = lat
    row, col = rs.index(x,y)
    value = rs.read(1)[row,col]
    return value

批量裁剪

import geopandas as gpd
import rasterio as rio
#from rasterio.mask import mask
import rasterio.mask  
from geopandas import GeoSeries
import os
def handle_grid(shpdatafile, rasterfile, out_file):
    """
    handle gird
    :param shpdatafile: path of shpfile
    :param rasterfile: path of rasterfile
    :param out_file
    :return:
    """
    shpdata = gpd.GeoDataFrame.from_file(shpdatafile)
    # print(shpdata.crs)
    # shpdata.crs = {'init': 'epsg:4326'}
    print(shpdata.crs)
    rasterdata = rio.open(rasterfile)
    print(rasterdata.crs)
    print(f"total {len(shpdata)}")
    for i in range(0, len(shpdata)):
        # getting vector data features
        print(f"current {i}th",end="")
        geo = shpdata.geometry[i]
        out_name = shpdata.gridid[i]
        #out_name = str(i)
        feature = [geo.__geo_interface__]
        # try:
        # Using feature to get building mask tif
        # print("running mask...")
        out_image, out_transform = rio.mask.mask(rasterdata, feature, all_touched=True, crop=True, nodata=rasterdata.nodata)
        #out_image, out_transform = rio.mask.mask(rasterdata, feature, all_touched=True, invert=True, nodata=rasterdata.nodata)
        # get tif Value值，and convert to list
        out_meta = rasterdata.meta.copy()
        out_meta.update({"driver": "GTiff",
                            "height": out_image.shape[1],
                            "width": out_image.shape[2],
                            "transform": out_transform})
        band_mask = rasterio.open(os.path.join(out_file, f"{out_name}.tif"), "w", **out_meta)
        # print(out_image)
        band_mask.write(out_image)
        # except Exception as e:
        #    print(e)
        #    pass
           
           
shp = "/home/data/jyx/population-5-5/run-data-result/urban_vitality/beijing_unicom/bj3857.shp"
raster = "/home/data/jyx/population-5-5/rs/uv_sentinel/beijing.tif"
out = "/home/data/jyx/population-5-5/rs/uv_sentinel/split_tif/"
handle_grid(shp,raster,out)