diff --git a/rewet/Output/Map.py b/rewet/Output/Map.py index b17e93b..bd158d4 100644 --- a/rewet/Output/Map.py +++ b/rewet/Output/Map.py @@ -22,8 +22,8 @@ import pandas as pd import numpy as np import shapely -import Output.Helper as Helper -import initial +import rewet.Output.Helper as Helper +#import initial #import time class Map(): @@ -33,34 +33,34 @@ def __init__(self): def loadShapeFile(self, shapeFileAddr='Northridge\GIS\Demand\demand_polygons.shp'): shape_file = gpd.read_file(shapeFileAddr) return shape_file - + def joinTwoShapeFiles(self, first, second): second = second.set_crs(crs=first.crs) joined_map = gpd.sjoin(first, second) - + return joined_map - + def createGeopandasPointDataFrameForNodes(self): s = gpd.GeoDataFrame(index=self.demand_node_name_list) point_list = [] point_name_list = [] - + for name in self.demand_node_name_list: coord = self.wn.get_node(name).coordinates point_list.append(shapely.geometry.Point(coord[0],coord[1])) point_name_list.append(name) s.geometry = point_list return s - + def getDLQNExceedenceProbabilityMap(self, data_frame, ihour , param): data = data_frame.transpose() scn_prob_list = self.scenario_prob #DLQN_dmg = pd.DataFrame(data=0, index=data.index, columns=data.columns) - + scn_prob = [scn_prob_list[scn_name] for scn_name in data.index] data['prob'] = scn_prob - + res_dict_list = [] tt = 0 if ihour: @@ -86,9 +86,9 @@ def getDLQNExceedenceProbabilityMap(self, data_frame, ihour , param): #min_ind = loop_dmg.idxmin() inter_value = loop_dmg.loc[min_ind, node_name] else: - + loop_dmg.loc['inter', 'ep'] = inter_ind - + loop_dmg = loop_dmg.sort_values('ep') ep_list = loop_dmg['ep'].to_list() inter_series = pd.Series(index=ep_list, data=loop_dmg[node_name].to_list()) @@ -102,7 +102,7 @@ def getDLQNExceedenceProbabilityMap(self, data_frame, ihour , param): else: for node_name in data_frame.index: loop_dmg = data[[node_name,'prob']] - + loop_dmg = loop_dmg.sort_values(node_name, ascending=False) loop_ep = Helper.EPHelper(loop_dmg['prob'].to_numpy()) loop_dmg['ep'] = loop_ep @@ -115,24 +115,24 @@ def getDLQNExceedenceProbabilityMap(self, data_frame, ihour , param): inter_value = loop_dmg.loc[min_ind, 'ep'] else: loop_dmg.loc['inter', node_name] = inter_ind - + loop_dmg = loop_dmg.sort_values(node_name) hour_list = loop_dmg[node_name].to_list() - + inter_series = pd.Series(index=hour_list, data=loop_dmg['ep'].to_list()) inter_series = inter_series.interpolate(method='linear') inter_value = inter_series.loc[inter_ind] if type(inter_value) != np.float64: inter_value = inter_value.mean() - + res_dict_list.append({'node_name':node_name, 'res':inter_value}) - + res = pd.DataFrame.from_dict(res_dict_list) res = res.set_index('node_name')['res'] - + s = self.createGeopandasPointDataFrameForNodes() s['res']=res - + #polygon = gpd.read_file('Northridge\GIS\Demand\demand_polygons.shp') #s = s.set_crs(epsg=polygon.crs.to_epsg()) #joined_map = gpd.sjoin(polygon, s) @@ -142,18 +142,18 @@ def getDLQNExceedenceProbabilityMap(self, data_frame, ihour , param): #props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) print(tt) return s - + def getOutageTimeGeoPandas_4(self, scn_name, LOS='DL' , iConsider_leak=False, leak_ratio=0, consistency_time_window=7200): #print(repr(LOS) + " " + repr(iConsider_leak)+" "+ repr(leak_ratio)+" "+repr(consistency_time_window ) ) self.loadScneariodata(scn_name) res = self.data[scn_name] map_res = pd.Series(data=0 , index=self.demand_node_name_list, dtype=np.int64) - + demands = self.getRequiredDemandForAllNodesandtime(scn_name) refined_res = res.node['demand'][self.demand_node_name_list] union_ = set(res.node['leak'].columns).union(set(self.demand_node_name_list) - (set(res.node['leak'].columns) ) - set(self.demand_node_name_list)) - (set(self.demand_node_name_list) - set(res.node['leak'].columns)) leak_res = res.node['leak'][union_] - + leak_data = [] if iConsider_leak: for name in leak_res: @@ -168,30 +168,30 @@ def getOutageTimeGeoPandas_4(self, scn_name, LOS='DL' , iConsider_leak=False, le leak_data = pd.DataFrame(leak_data).transpose() demands = demands[self.demand_node_name_list] - + if LOS=="DL": DL_res_not_met_bool = refined_res <= demands * 0.01 elif LOS =="QN": DL_res_not_met_bool = refined_res < demands * 0.98 - + time_window = consistency_time_window + 1 time_list = DL_res_not_met_bool.index.to_list() time_list.reverse() - + for time in time_list: past_time_beg = time - time_window window_data = DL_res_not_met_bool.loc[past_time_beg:time] window_data = window_data.all() window_data_false = window_data[window_data == False] DL_res_not_met_bool.loc[time, window_data_false.index] = False - + for name in DL_res_not_met_bool: if name in leak_data.columns: leak_data_name = leak_data[name] for time in leak_data_name.index: if leak_data_name.loc[time] == True: DL_res_not_met_bool.loc[time, name] = True - + all_node_name_list = refined_res.columns only_not_met_bool = DL_res_not_met_bool.any(0) only_not_met_any = all_node_name_list[only_not_met_bool] @@ -203,7 +203,7 @@ def getOutageTimeGeoPandas_4(self, scn_name, LOS='DL' , iConsider_leak=False, le rolled_DL_res_MET = DL_res_MET[name].rolling(time_window, center=True).sum() rolled_DL_res_MET = rolled_DL_res_MET.sort_index(ascending=False) rolled_DL_res_MET.dropna(inplace=True) - + false_found, found_index = Helper.helper_outageMap(rolled_DL_res_MET.ge(time_window-1)) #if name == "SM323": #return DL_res_MET[name], rolled_DL_res_MET, false_found, rolled_DL_res_MET.index[found_index], rolled_DL_res_MET.ge(time_window-1), found_index @@ -215,42 +215,42 @@ def getOutageTimeGeoPandas_4(self, scn_name, LOS='DL' , iConsider_leak=False, le else: latest_time = rolled_DL_res_MET.index[found_index] latest_time = rolled_DL_res_MET.index[found_index] - + map_res.loc[name] = latest_time #map_res = map_res/(3600*24) geopandas_df = self.createGeopandasPointDataFrameForNodes() geopandas_df.loc[map_res.index.to_list(), 'restoration_time'] = map_res.to_list() - + return geopandas_df - + def getOutageTimeGeoPandas_5(self, scn_name, bsc='DL' , iConsider_leak=False, leak_ratio=0, consistency_time_window=7200, sum_time=False): self.loadScneariodata(scn_name) res = self.data[scn_name] map_res = pd.Series(data=0 , index=self.demand_node_name_list, dtype=np.int64) - + required_demand = self.getRequiredDemandForAllNodesandtime(scn_name) delivered_demand = res.node['demand'][self.demand_node_name_list] common_nodes_leak = list (set( res.node['leak'].columns ).intersection( set( self.demand_node_name_list ) )) leak_res = res.node['leak'][common_nodes_leak] - + common_nodes_demand = list( set(delivered_demand.columns).intersection(set(self.demand_node_name_list) ) ) delivered_demand = delivered_demand[common_nodes_demand] required_demand = required_demand[common_nodes_demand] - + required_demand.sort_index(inplace=True) delivered_demand.sort_index(inplace=True) leak_res.sort_index(inplace=True) - + #return delivered_demand, required_demand, leak_res - + if bsc=="DL": bsc_res_not_met_bool = (delivered_demand.fillna(0) <= required_demand * 0.0001) elif bsc =="QN": bsc_res_not_met_bool = (delivered_demand.fillna(0) < required_demand * 0.9999) else: raise ValueError("Unknown BSC= "+str(bsc)) - + if iConsider_leak : #return leak_res, required_demand leak_res_non_available_time_list = set(required_demand[leak_res.columns].index) - set(leak_res.index) @@ -264,15 +264,15 @@ def getOutageTimeGeoPandas_5(self, scn_name, bsc='DL' , iConsider_leak=False, le combined_negative_result = (bsc_res_not_met_bool | leak_criteria_exceeded).dropna(axis=1) #return combined_negative_result bsc_res_not_met_bool.loc[:, combined_negative_result.columns] = combined_negative_result - + #end_time = delivered_demand.min() end_time = delivered_demand.index.max() if consistency_time_window > 1: time_beg_step_list = np.arange(0, end_time, consistency_time_window) - + #time_beg_step_list = np.append(time_beg_step_list, [end_time]) time_end_step_list = time_beg_step_list # + consistency_time_window - window_bsc_not_met = pd.DataFrame(index=time_end_step_list, columns= bsc_res_not_met_bool.columns, dtype=bool) + window_bsc_not_met = pd.DataFrame(index=time_end_step_list, columns= bsc_res_not_met_bool.columns, dtype=bool) #return bsc_res_not_met_bool#, delivered_demand, required_demand for step_time_beg in time_beg_step_list: step_time_end = step_time_beg + consistency_time_window @@ -285,19 +285,19 @@ def getOutageTimeGeoPandas_5(self, scn_name, bsc='DL' , iConsider_leak=False, le window_bsc_not_met.drop(step_time_beg, inplace=True) else: window_bsc_not_met = bsc_res_not_met_bool - + pre_incident = (window_bsc_not_met.loc[:3600*3]).any() non_incident = pre_incident[pre_incident==False].index - + not_met_node_name_list = window_bsc_not_met.any() - + #("****************") #print(not_met_node_name_list[not_met_node_name_list==True]) - + not_met_node_name_list = not_met_node_name_list[not_met_node_name_list==True] not_met_node_name_list = not_met_node_name_list.index - - + + if sum_time: time_difference = window_bsc_not_met.index[1:] - window_bsc_not_met.index[:-1] timed_diference_window_bsc_not_met = \ @@ -305,48 +305,48 @@ def getOutageTimeGeoPandas_5(self, scn_name, bsc='DL' , iConsider_leak=False, le timed_diference_window_bsc_not_met.iloc[0] = 0 sum_window_bsc_not_met = timed_diference_window_bsc_not_met.sum() return sum_window_bsc_not_met - + window_bsc_not_met = window_bsc_not_met[not_met_node_name_list] cut_time = window_bsc_not_met.index.max() non_incident = list( set(non_incident).intersection(set(not_met_node_name_list) ) ) for step_time, row in window_bsc_not_met[non_incident].iterrows(): if step_time <= 14400: continue - + if row.any() == False: print(step_time) cut_time = step_time break - + window_bsc_not_met = window_bsc_not_met.loc[:cut_time] window_bsc_not_met = window_bsc_not_met.loc[:cut_time] - + #return window_bsc_not_met #print(not_met_node_name_list) time_bsc_not_met_time = window_bsc_not_met.sort_index(ascending=False).idxmax() map_res.loc[time_bsc_not_met_time.index] = time_bsc_not_met_time - - + + never_reported_nodes = set(self.demand_node_name_list) - set(common_nodes_demand) number_of_unreported_demand_nodes = len(never_reported_nodes) if number_of_unreported_demand_nodes > 0: warnings.warn("REWET WARNING: there are " + str(number_of_unreported_demand_nodes ) + "unreported nodes") map_res.loc[never_reported_nodes ] = end_time - + map_res = map_res/(3600*24) return map_res - + s = gpd.GeoDataFrame(index=self.demand_node_name_list) point_list = [] point_name_list = [] - + for name in self.demand_node_name_list: coord = self.wn.get_node(name).coordinates point_list.append(shapely.geometry.Point(coord[0],coord[1])) point_name_list.append(name) - + s.geometry = point_list - + s.loc[map_res.index.to_list(), 'restoration_time'] = map_res.to_list() polygon = gpd.read_file('Northridge\GIS\Demand\demand_polygons.shp') @@ -354,36 +354,36 @@ def getOutageTimeGeoPandas_5(self, scn_name, bsc='DL' , iConsider_leak=False, le joined_map = gpd.sjoin(polygon, s) #return joined_map #joined_map.loc[map_res.index.to_list(), 'restoration_time'] = (map_res/3600/24).to_list() - + return joined_map - + def percentOfEffectNodes(self, scn_name, bsc='QN' , iConsider_leak=True, leak_ratio=0.75, consistency_time_window=7200): self.loadScneariodata(scn_name) res = self.data[scn_name] map_res = pd.Series(data=0 , index=self.demand_node_name_list, dtype=np.int64) - + required_demand = self.getRequiredDemandForAllNodesandtime(scn_name) delivered_demand = res.node['demand'][self.demand_node_name_list] common_nodes_leak = set(res.node['leak'].columns).intersection(set(self.demand_node_name_list)) leak_res = res.node['leak'][common_nodes_leak] - + common_nodes_demand = list( set(delivered_demand.columns).intersection(set(self.demand_node_name_list) ) ) delivered_demand = delivered_demand[common_nodes_demand] required_demand = required_demand[common_nodes_demand] - + required_demand.sort_index(inplace=True) delivered_demand.sort_index(inplace=True) leak_res.sort_index(inplace=True) - + #return delivered_demand, required_demand, leak_res - + if bsc=="DL": bsc_res_not_met_bool = (delivered_demand.fillna(0) <= required_demand * 0.1) elif bsc =="QN": bsc_res_not_met_bool = (delivered_demand.fillna(0) < required_demand * 0.98) else: raise ValueError("Unknown BSC= "+str(bsc)) - + if iConsider_leak : #return leak_res, required_demand leak_res_non_available_time_list = set(required_demand[leak_res.columns].index) - set(leak_res.index) @@ -397,14 +397,14 @@ def percentOfEffectNodes(self, scn_name, bsc='QN' , iConsider_leak=True, leak_ra combined_negative_result = (bsc_res_not_met_bool | leak_criteria_exceeded).dropna(axis=1) #return combined_negative_result bsc_res_not_met_bool.loc[:, combined_negative_result.columns] = combined_negative_result - + #end_time = delivered_demand.min() end_time = delivered_demand.index.max() time_beg_step_list = np.arange(0, end_time, consistency_time_window) - + #time_beg_step_list = np.append(time_beg_step_list, [end_time]) time_end_step_list = time_beg_step_list # + consistency_time_window - window_bsc_not_met = pd.DataFrame(index=time_end_step_list, columns= bsc_res_not_met_bool.columns, dtype=bool) + window_bsc_not_met = pd.DataFrame(index=time_end_step_list, columns= bsc_res_not_met_bool.columns, dtype=bool) #return bsc_res_not_met_bool#, delivered_demand, required_demand for step_time_beg in time_beg_step_list: step_time_end = step_time_beg + consistency_time_window @@ -418,48 +418,48 @@ def percentOfEffectNodes(self, scn_name, bsc='QN' , iConsider_leak=True, leak_ra #return window_bsc_not_met pre_incident = (window_bsc_not_met.loc[:3600*3]).any() non_incident = pre_incident[pre_incident==False].index - + number_of_good_nodes = len(non_incident) - + not_met_node_name_list = window_bsc_not_met.any() - + #("****************") #print(not_met_node_name_list[not_met_node_name_list==True]) - + not_met_node_name_list = not_met_node_name_list[not_met_node_name_list==True] not_met_node_name_list = not_met_node_name_list.index window_bsc_not_met = window_bsc_not_met[not_met_node_name_list] - + cut_time = window_bsc_not_met.index.max() non_incident = list( set(non_incident).intersection(set(not_met_node_name_list) ) ) for step_time, row in window_bsc_not_met[non_incident].iterrows(): if step_time <= 14400: continue - + if row.any() == False: print(step_time) cut_time = step_time break - + cut_time = 24 * 3600 window_bsc_not_met = window_bsc_not_met.loc[:cut_time] window_bsc_not_met = window_bsc_not_met.loc[:cut_time] - - + + number_of_bad_node_at_damage = window_bsc_not_met[non_incident].loc[14400].sum() percent_init = number_of_bad_node_at_damage / number_of_good_nodes * 100 #return window_bsc_not_met #print(not_met_node_name_list) time_bsc_not_met_time = window_bsc_not_met.sort_index(ascending=False).idxmax() map_res.loc[time_bsc_not_met_time.index] = time_bsc_not_met_time - - + + never_reported_nodes = set(self.demand_node_name_list) - set(common_nodes_demand) number_of_unreported_demand_nodes = len(never_reported_nodes) if number_of_unreported_demand_nodes > 0: warnings.warn("REWET WARNING: there are " + str(number_of_unreported_demand_nodes ) + "unreported nodes") map_res.loc[never_reported_nodes ] = end_time - + map_res = map_res/(3600*24) percent = (map_res.loc[non_incident] > 0).sum() / number_of_good_nodes * 100 - return np.round(percent_init, 2) , np.round(percent, 2) \ No newline at end of file + return np.round(percent_init, 2) , np.round(percent, 2)