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: from geopandas.tools import overlay data = 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()