153 lines
4.2 KiB
Python
153 lines
4.2 KiB
Python
import numpy as np
|
|
import os
|
|
from geopandas import *
|
|
from pandas import *
|
|
from shapely.geometry import Point
|
|
import rasterio as rio
|
|
import rasterio.mask
|
|
from ..util.utils import *
|
|
|
|
|
|
def get_data(path):
|
|
data = []
|
|
with open(path, 'r') as file:
|
|
for line in file.readlines():
|
|
arr = line.split()
|
|
data.append(float(arr[4]))
|
|
return np.asarray(data)
|
|
|
|
|
|
def get_stations():
|
|
data = []
|
|
for i in range(1, 180):
|
|
path = r"../data/rain/%s.txt" % str(i).rjust(3, '0')
|
|
with open(path, 'r') as file:
|
|
first_line = file.readline(100)
|
|
arr = first_line.split()
|
|
data.append(str(i).rjust(3, '0') + ',' + arr[0] + ',' + arr[1] + '\r')
|
|
with open('../data/rain/station.csv', 'w') as file:
|
|
file.writelines(data)
|
|
|
|
|
|
def get_station_polygon():
|
|
path = '../data/station/station_new.csv'
|
|
data = pandas.read_csv(path)
|
|
geometry = [Point(xy) for xy in zip(data.X, data.Y)]
|
|
data['geometry'] = geometry
|
|
return get_buffer_data(GeoDataFrame(data), dis=10000)
|
|
|
|
|
|
def read_from_out_path(type, out_path):
|
|
if type == 'spi':
|
|
index = 3
|
|
else:
|
|
index = 4
|
|
with open(out_path, "r") as out_file:
|
|
lines = out_file.readlines()[index:]
|
|
|
|
def convert(str):
|
|
try:
|
|
return float(str.replace('\n', ''))
|
|
except:
|
|
return 1
|
|
|
|
lines = [convert(line) for line in lines]
|
|
data = np.asarray(lines)
|
|
return data
|
|
|
|
|
|
def prepare_data_file(path, input_path, type, forecast_data=None):
|
|
with open(path, 'r') as file1:
|
|
with open(input_path, 'w') as file2:
|
|
file2_lines = []
|
|
file1_lines = file1.readlines()
|
|
arr = file1_lines[0].split(',')
|
|
year, month = arr[2], arr[3]
|
|
file2_lines.append(os.path.basename(path)[: os.path.basename(path).index('.')] + "\n")
|
|
if type == 'spei':
|
|
file2_lines.append(arr[1] + "\n")
|
|
file2_lines.append(year + ";" + month + "\n")
|
|
file2_lines.append('12\n')
|
|
for line in file1_lines:
|
|
arr = line.split(',')
|
|
if type == 'spei':
|
|
file2_lines.append(arr[4] + ';' + arr[5].replace('\r', '\n'))
|
|
else:
|
|
file2_lines.append(arr[4] + '\n')
|
|
if forecast_data:
|
|
forecast_data = [str(d) + '\n' for d in forecast_data]
|
|
file2_lines.extend(forecast_data)
|
|
file2.writelines(file2_lines)
|
|
|
|
|
|
def forecast_by_exe(path, scale, forecast_data, type='spi'):
|
|
spi_data = calc_spi_by_exe(path, scale, type, forecast_data)
|
|
return spi_data[-1 * len(forecast_data):]
|
|
|
|
|
|
def calc_spi_by_exe(path, scale, type='spi', forecast_data=None):
|
|
if type == 'spi':
|
|
spi_path = 'spi.exe'
|
|
else:
|
|
spi_path = 'spei.exe'
|
|
input_path = "../data/temp/input.csv"
|
|
prepare_data_file(path, input_path, type, forecast_data)
|
|
|
|
out_path = "../data/temp/res.csv"
|
|
args = [spi_path, str(scale), input_path, out_path]
|
|
call_exe(args)
|
|
return read_from_out_path(type, out_path)
|
|
|
|
|
|
def call_exe(args):
|
|
import subprocess
|
|
p1 = subprocess.Popen(args, stdout=subprocess.PIPE)
|
|
# print(p1.communicate()[0])
|
|
p1.wait()
|
|
|
|
|
|
def get_spi_dataframe(station_polygon, newdata_path):
|
|
spi_list = []
|
|
spei_list = []
|
|
|
|
for station in range(1, 334):
|
|
path = r"../data/station/station_new/V7%s.csv" % str(station).rjust(3, '0')
|
|
if not os.path.exists(path):
|
|
continue
|
|
|
|
forecast_spi = forecast_by_exe(path, 1, ["50", "3", "12"], 'spi') # spi
|
|
forecast_spei = forecast_by_exe(path, 1, ["2;-2", "10;3", "10;5"], 'spei') # spei
|
|
spi_list.append(forecast_spi[-1])
|
|
spei_list.append(forecast_spei[-1])
|
|
station_polygon['spi'] = spi_list
|
|
station_polygon['spei'] = spei_list
|
|
return station_polygon
|
|
|
|
|
|
def get_disaster_gdf(gdf, type):
|
|
return gdf[gdf[type] < 0]
|
|
|
|
|
|
if __name__ == '__main__':
|
|
station_buffer = get_station_polygon()
|
|
spi_gdf = get_spi_dataframe(station_buffer, '')
|
|
# SPI
|
|
spi_disaster = get_disaster_gdf(spi_gdf, 'spi')
|
|
# spi_disaster.to_file('spi.shp')
|
|
res_spi = agg_raster_by_gdf_value('../data/raster/gdp_3857.tif', spi_disaster)
|
|
print(res_spi)
|
|
# SPI 三区九县
|
|
county_spi = intersection_xzqh(spi_disaster)
|
|
res_county_spi = agg_raster_by_gdf_value('../data/raster/gdp_3857.tif', county_spi)
|
|
print(res_county_spi)
|
|
|
|
# SPEI
|
|
spei_disaster = get_disaster_gdf(spi_gdf, 'spei')
|
|
res_spei = agg_raster_by_gdf_value('../data/raster/gdp_3857.tif', spei_disaster)
|
|
print(res_spei)
|
|
|
|
# SPEI 三区九县
|
|
county_spei = intersection_xzqh(spei_disaster)
|
|
res_county_spei = agg_raster_by_gdf_value('../data/raster/gdp_3857.tif', county_spei)
|
|
print(res_county_spei)
|