drought/disaster/util/utils.py

from geopandas import *
from pandas import *
from shapely.geometry import Point
import rasterio as rio
import rasterio.mask
from rasterio.warp import reproject, Resampling

DATA_BASE_PATH = '/home/g214/data_from_chenhao/data_analyse/'
# DATA_BASE_PATH = r"D:\7. business\7.baoji\code\drought_analyse\data/"

def get_buffer_data(gdf, crs={'init': 'epsg:4326', 'no_defs': True}, dis=1):
	data = gdf
	if not data.crs:
		data.crs = crs
	else:
		crs = data.crs

	data = data.to_crs(epsg=3857)
	buffer = data.buffer(dis)
	data.set_geometry(buffer, inplace=True)
	data.to_crs(crs=crs, inplace=True)
	return data


def read_xzqh(path):
	shp_bjqh = GeoDataFrame.from_file(path)
	shp_bjqh = shp_bjqh[['geometry', 'county']]
	# shp_bjqh.rename(columns={'县名称': 'county'}, inplace=True)
	shp_bjqh.to_crs(epsg=4326, inplace=True)
	return shp_bjqh


def intersection(gdf1, gdf2):
	try:
		data = geopandas.overlay(gdf1, gdf2, how='intersection')
		data.crs = gdf1.crs
		return data
	except Exception as e:
		print(str(e))
		return None


def intersection_xzqh(gdf):
	bjqh = read_xzqh(DATA_BASE_PATH + 'boundary/bj.shp')
	intersection_data = intersection(gdf, bjqh)
	return intersection_data


def join(gdf1, gdf2):
	return geopandas.sjoin(gdf1, gdf2, how="inner", op='intersects')


def agg_raster_by_gdf_value(path, gdf):
	# crs = '+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
	with rio.open(path) as src:
		gdf = gdf.to_crs(epsg=3857)
		sum = {}
		for i in range(0, len(gdf)):
			geo = gdf.iloc[i].geometry
			try:
				feature = [GeoSeries(geo).__geo_interface__['features'][0]['geometry']]
				out_image, out_transform = rio.mask.mask(src, feature, crop=True, nodata=src.nodata)
				res = np.nansum(np.where(out_image == src.nodata, np.nan, out_image))
				if 'county' in gdf:
					county = gdf.iloc[i].county
					print(county)
					if county in sum.keys():
						sum[county] = sum.get(county) + res
					else:
						sum[county] = res
				elif 'all' in sum.keys():
					sum['all'] = sum.get('all') + res
				else:
					sum['all'] = res
			# for band in out_image:
			#     for row in band:
			#         for column in row:
			#             if column != src.nodata:
			#                 sum += column
			# print(sum)
			except Exception as e:
				print(str(e))
		return sum


def agg_raster_by_gdf_area(path, gdf):
	# crs = '+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
	with rio.open(path) as src:
		gdf = gdf.to_crs(epsg=3857)
		sum = {}
		for i in range(0, len(gdf)):
			geo = gdf.iloc[i].geometry
			try:
				feature = [GeoSeries(geo).__geo_interface__['features'][0]['geometry']]
				out_image, out_transform = rio.mask.mask(src, feature, crop=True, nodata=src.nodata)
				out_image = np.where(out_image == 1, 1, np.nan)
				print(out_image)
				res = np.nansum(out_image) * src.width * src.height
				if 'county' in gdf:
					county = gdf.iloc[i].county
					if county in sum.keys():
						sum[county] = sum.get(county) + res
					else:
						sum[county] = res
				elif 'all' in sum.keys():
					sum['all'] = sum.get('all') + res
				else:
					sum['all'] = res
			except:
				pass
		return sum


def agg_prevent_disater_by_gdf_area(path, gdf):
	# crs = '+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
	with rio.open(path) as src:
		gdf = gdf.to_crs(epsg=3857)
		sum = {}
		for i in range(0, len(gdf)):
			geo = gdf.iloc[i].geometry
			try:
				feature = [GeoSeries(geo).__geo_interface__['features'][0]['geometry']]
				out_image, out_transform = rio.mask.mask(src, feature, crop=True, nodata=src.nodata)
				out_image = np.where(out_image == 1, 1, np.nan)
				print(out_image)
				res = np.nansum(out_image) * src.width * src.height
				name = str(gdf.iloc[i].name)
				if name in sum.keys():
					sum[name] = sum.get(name) + res
				else:
					sum[name] = res
			except:
				pass
		return sum


def agg_gdf_by_gdf_area(gdf1, gdf2):
	intersection_data = intersection(gdf1, gdf2)
	area = intersection_data.to_crs(epsg=3857).area
	intersection_data['area'] = area
	intersection_data = intersection_data[['name', 'area']]
	# print(intersection_data.head())
	return intersection_data.groupby('name')['area'].sum()


def agg_raster_by_gdf_area_and_aspect_slope(path, path_aspect, path_slope, gdf):
	with rio.open(path) as src:
		with rio.open(path_aspect) as src1:
			with rio.open(path_slope) as src2:
				gdf = gdf.to_crs(epsg=3857)
				res_flat = 0
				res_30 = 0
				res_more_than_30 = 0
				for i in range(0, len(gdf)):
					geo = gdf.iloc[i].geometry
					try:
						feature = [GeoSeries(geo).__geo_interface__['features'][0]['geometry']]
						out_image, out_transform = rio.mask.mask(src, feature, crop=True, nodata=src.nodata)
						out_image_aspect, out_transform_aspect = rio.mask.mask(src1, feature, crop=True,
						                                                       nodata=src1.nodata)
						out_image_slope, _ = rio.mask.mask(src2, feature, crop=True, nodata=src2.nodata)
						out_image_aspect_reproject = np.empty(shape=(out_image.shape[0],  # same number of bands
						                                 round(out_image.shape[1]),
						                                 round(out_image.shape[2])))

						out_image_slope_reproject = np.empty(shape=(out_image.shape[0],  # same number of bands
						                                 round(out_image.shape[1]),
						                                 round(out_image.shape[2])))
						reproject(
							out_image_aspect, out_image_aspect_reproject,
							src_transform=out_transform_aspect,
							dst_transform=out_transform,
							src_crs='+proj=latlong',
							dst_crs='+proj=latlong',
							resampling=Resampling.bilinear)
						reproject(
							out_image_slope, out_image_slope_reproject,
							src_transform=out_transform_aspect,
							dst_transform=out_transform,
							# src_crs={'init': 'EPSG:3857'},
							# dst_crs={'init': 'EPSG:3857'},
							resampling=Resampling.bilinear)
						out_image_slope_reproject = np.where(out_image_slope_reproject == src1.nodata, 1000, out_image_slope_reproject)
						out_image_slope_reproject = out_image_slope_reproject.astype(int)
						out_image_aspect_reproject = np.where(out_image_aspect_reproject == src2.nodata, 1000, out_image_aspect_reproject)
						out_image_aspect_reproject = out_image_aspect_reproject.astype(int)
						for i in range(0, out_image.shape[1]):
							for j in range(0, out_image.shape[2]):
								aspect = out_image_aspect_reproject[0, i, j]
								slope = out_image_slope_reproject[0, i, j]
								print(aspect, slope, out_image[0, i, j])
								if out_image[0, i, j] > 0:
									if slope <= 0:
										print('slope')
										if aspect < 0:
											res_flat += 1
									elif slope < 30:
										res_30 += 1
									else:
										res_more_than_30 += 1
						print(res_flat, res_30, res_more_than_30)
					except Exception as e:
						print(str(e))
				return (res_flat, res_30, res_more_than_30)

def merge_gdf_geo_to_geojson(gdf):
	return gdf.to_json()