把pandas的dataframe转为geopandas的地理格式（df to geodf）
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
重采样
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)
log_resolver.ipynb用来绘图
原因是通过程序绘制的图更规范，有一些指标，如散点图的拟合公式、\(R^2\)等 ，以前在excel中需要点点点才能显示，在观察、分析大量log文件时，很繁琐，不利于发现规律。
而在程序写好之后，log可以很直观的看到，并进行对比。并且由于程序制图较规范，便于不同实验之间进行比较，美观度也还可以，后续略作调整即可用于论文插图
以下代码实现了plt的density scatter.
def density_scatter( x , y, ax = None, sort = True, bins = 20, is_cbar=True, **kwargs ) :
"""
Scatter plot colored by 2d histogram
"""
if ax is None :
fig , ax = plt.subplots()
data , x_e, y_e = np.histogram2d( x, y, bins = bins, density = True )
z = interpn( ( 0.5*(x_e[1:] + x_e[:-1]) , 0.5*(y_e[1:]+y_e[:-1]) ) , data , np.vstack([x,y]).T , method = "splinef2d", bounds_error = False)
#To be sure to plot all data
z[np.where(np.isnan(z))] = 0.0
# Sort the points by density, so that the densest points are plotted last
if sort :
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
ax.scatter( x, y, c=z, **kwargs )
ax.grid(linestyle='--',linewidth=0.5)
norm = Normalize(vmin = np.min(z), vmax = np.max(z))
if is_cbar:
cbar = plt.colorbar(cm.ScalarMappable(norm = norm), ax=ax)
cbar.ax.set_ylabel('Density')
return ax