{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Deriving all KPIs related to cold/winter weather for cities"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nfrom scipy import stats\nfrom sklearn.preprocessing import RobustScaler, QuantileTransformer\nimport rasterio\nfrom matplotlib import colors\nfrom scipy import stats\nfrom tqdm import *\nfrom pyproj import Transformer\nimport ast\nimport glob\nfrom io import StringIO","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"path = \"../input/externald/Supplementary Data/\" \nvulnerability = pd.read_csv(path + \"cities_updated_geo_us_2020.csv\")\nacct_number = vulnerability[\"Account Number\"]\ncounties = (vulnerability['counties'])\ncounties_weights = (vulnerability['counties_weights'])\nstates = (vulnerability['state'])\nflag_multi_counties = vulnerability['flag_multiple_counties']\nmatch_fips = vulnerability['FIPS']\nMSA = (vulnerability['MSA'])\nmatch_fips = [ast.literal_eval(match_fips[i]) if flag_multi_counties[i] else int(match_fips[i]) for i in range(len(match_fips))]\ncities = vulnerability[\"city_bing\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cities_2020 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2020_Full_Cities_Dataset.csv\")\n\ncities_disc_2020 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2020_Cities_Disclosing_to_CDP.csv\")\n\ncities_2020_merged = pd.merge(cities_2020, cities_disc_2020, on=\"Account Number\", how=\"outer\", copy=False)\nunique_20 = np.unique(cities_disc_2020[cities_disc_2020[\"Country\"] == \"United States of America\"][\"Account Number\"])\ncities_2020_merged_sub = cities_2020_merged[np.in1d(cities_2020_merged[\"Account Number\"], unique_20)]\n\ndef fetch_answer_individual(account, question_number, org_type=\"city\", year=2018, corp_res_type=None, column_number=None, row_number=None):\n    if(org_type==\"city\"):\n        df = cities_2020_merged_sub\n        subset = df[df['Question Number'] == question_number]\n        answer = subset[subset[\"Account Number\"] == account]\n        if((column_number is not None) and (row_number is None)):\n            answer = answer[answer[\"Column Number\"] == column_number]\n        elif((column_number is not None) and (row_number is not None)):\n            answer = answer[(answer[\"Column Number\"] == column_number) & (answer[\"Row Number\"] == row_number)]\n        elif((column_number is None) and (row_number is not None)):\n            answer = answer[answer[\"Row Number\"] == row_number]\n        else:\n            pass\n    elif(org_type==\"corp\" and corp_res_type==\"cc\"):\n        df = corporations_2020\n        subset = df[df['question_number'] == question_number]\n        answer = subset[subset[\"account_number\"] == account]\n        if((column_number is not None) and (row_number is None)):\n            answer = answer[answer[\"column_number\"] == column_number]\n        elif((column_number is not None) and (row_number is not None)):\n            answer = answer[(answer[\"column_number\"] == column_number) & (answer[\"row_number\"] == row_number)]\n        elif((column_number is None) and (row_number is not None)):\n            answer = answer[answer[\"row_number\"] == row_number]\n        else:\n            pass\n    elif(org_type==\"corp\" and corp_res_type==\"ws\"):\n        df = all_corps_cc[year]\n        subset = df[df['question_number'] == question_number]\n        answer = subset[subset[\"account_number\"] == account]\n    else:\n        print(\"Something went wrong. Try again.\")\n    \n    return answer","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score_hazard_risk(cdf):\n    prob = scores_lo_hi[str(cdf[cdf[\"Column Number\"] == 3.0][\"Response Answer\"].iloc[0])]\n    mag = scores_lo_hi[str(cdf[cdf[\"Column Number\"] == 4.0][\"Response Answer\"].iloc[0])]\n    \n    future_freq = scores_increasing_decreasing[str(cdf[cdf[\"Column Number\"] == 8.0][\"Response Answer\"].iloc[0])]\n    future_intensity = scores_increasing_decreasing[str(cdf[cdf[\"Column Number\"] == 9.0][\"Response Answer\"].iloc[0])]\n    \n    future_mag = scores_lo_hi[str(cdf[cdf[\"Column Number\"] == 10.0][\"Response Answer\"].iloc[0])]\n    return (prob*mag) + (prob*future_freq*future_intensity*future_mag)\n\ndef return_risk_and_score_per_county(hazard_type):\n    risk = []\n    self_reported_score = []\n\n    for i in range(len(vulnerability)):\n        df = fetch_answer_individual(vulnerability[\"Account Number\"].iloc[i], '2.1')\n        idx = np.where(np.array(df[\"Response Answer\"]) == hazard_type)[0]\n        if(len(idx) > 0):\n            risk.append(1)\n            cdf = df[df[\"Row Number\"] == df[\"Row Number\"].iloc[idx[0]]]\n            self_reported_score.append(score_hazard_risk(cdf))\n        else:\n            risk.append(np.nan)\n            self_reported_score.append(np.nan)\n    return risk, np.nan_to_num(np.array(self_reported_score))\n\nscores_lo_hi = {\n    \"Low\" : 1, \n    \"Medium Low\": 2, \n    \"Medium\":3,\n    \"Medium High\":4,\n    \"High\":5,\n    \"Do not know\":0.5,\n    \"Does not currently impact the city\":0.5,\n    \"nan\":0,\n}\n\nscores_increasing_decreasing = {\n    \"Increasing\" : 2,\n    \"Decreasing\" : 0.5,\n    \"None\" : 1, \n    \"Do not know\" : 1, \n    \"Not expected to happen in the future\" : 0.1,\n    \"nan\" : 0.0\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Scoring survey response \n\nWe can derive survey score using the equation:\n#### (probability x magnitude) + (probability x future_frequency x future_intensity x future_mag)"},{"metadata":{"trusted":true},"cell_type":"code","source":"risk_cold_wave, score_cold_wave = return_risk_and_score_per_county(\"Extreme cold temperature > Cold wave\")\nrisk_cold_days, score_cold_days = return_risk_and_score_per_county(\"Extreme cold temperature > Extreme cold days\")\nrisk_winter_conditions, score_winter_conditions = return_risk_and_score_per_county(\"Extreme cold temperature > Extreme winter conditions\")\nrisk_heavy_snow, score_heavy_snow = return_risk_and_score_per_county(\"Extreme Precipitation > Heavy snow\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def standardize_rank(arr, direction=1):\n    rank = (stats.mstats.rankdata(np.ma.masked_invalid(arr)))\n    rank[rank == 0] = np.nan\n    if(direction == 1):\n        return rank/(np.nanmax(rank)) \n    else:\n        return 1-(rank/(np.nanmax(rank)))\n\n    \ncold_wave_ranked = standardize_rank(score_cold_wave, direction=1)\ncold_days_ranked = standardize_rank(score_cold_wave, direction=1)\nwinter_conditions_ranked = standardize_rank(score_cold_wave, direction=1)\nheavy_snow_ranked = standardize_rank(score_cold_wave, direction=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fema = pd.read_csv(\"../input/femadisasters/DisasterDeclarationsSummaries.csv\")\n\nct = np.char.zfill(np.array(fema['fipsCountyCode']).astype(str), 3)\nst = np.char.zfill(np.array(fema['fipsStateCode']).astype(str), 2)\nfema['fips'] = np.core.defchararray.add(st, ct).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fema_cold_events = []\nfor i in range(len(match_fips)):\n    if(not flag_multi_counties[i]):\n        fips = int(match_fips[i])\n        cdf = fema[(fema[\"incidentType\"] == \"Severe Ice Storm\") | (fema[\"incidentType\"] == \"Snow\") | (fema[\"incidentType\"] == \"Freezing\")]\n        cdf = cdf[cdf[\"fips\"] == fips]\n        fema_cold_events.append(len(cdf))\n    else:\n        fips = match_fips[i]\n        temp = []\n        for j in range(len(fips)):\n            cdf = fema[(fema[\"incidentType\"] == \"Severe Ice Storm\") | (fema[\"incidentType\"] == \"Snow\")| (fema[\"incidentType\"] == \"Freezing\")]\n            cdf = cdf[cdf[\"fips\"] == fips[j]]\n            temp.append(len(cdf))\n        ar = np.array(temp)\n        weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n        fema_cold_events.append(weighted_avg)\n\nfema_cold_events_ranked = standardize_rank(fema_cold_events, direction=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"str_tmin = open('../input/tmindatasetnoaa/climdiv-tmincy-v1.0.0-20201104', 'r').read()\ncolumns = [\"Code\", \"Jan\", \"Feb\", \"Mar\", \"Apr\", \"May\", \"Jun\", \"Jul\", \"Aug\", \"Sep\", \"Oct\", \"Nov\", \"Dec\"]\nfor df in pd.read_fwf(StringIO(str_tmin), header=None, chunksize=500000, names=columns,\n                      converters={h:str for h in columns}):\n    tmin = df\n\nyears = []\nfips = []\nfor i in range(len(tmin)):\n    fips.append(int(str(tmin[\"Code\"].iloc[i])[0:5]))\n    years.append(int(str(tmin[\"Code\"].iloc[i])[7:11]))\n    \ntmin[\"year\"] = np.array(years).astype(int)\ntmin[\"fips\"] = np.array(fips).astype(int)\nstate_fips = pd.read_csv(\"../input/statefips/state_fips.csv\")\ndf_states = state_fips.drop([1,10])\ndf_states = df_states.reset_index()\ndf_states[\"index\"] = np.array(df_states.index) + 1\nfips_to_code = dict(zip(df_states[\"fips\"], df_states[\"index\"]))\nfips_to_code[15] = 99 # HI\nfips_to_code[11] = 99 # DC\nfips_to_code[2] = 99 # AK","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tmin_averaged = []\nfor i in range(len(match_fips)):\n    if(not flag_multi_counties[i]):\n        fips = int(match_fips[i])\n        padstr = str(np.char.zfill(str(fips), width=5))\n        st_code = int((padstr)[0:2])\n        fips = int(str(fips_to_code[st_code]) + padstr[2:])\n        cdf = tmin[(tmin['fips'] == fips) & (tmin['year'] > 2000)]\n        cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float)  \n        cdf = cdf[(cdf != 0.00) & (cdf != -9.99)]\n        tmin_averaged.append(np.mean(np.min(cdf, axis=1)))\n    else:\n        fips = match_fips[i]\n        temp = []\n        for j in range(len(fips)):\n            padstr = str(np.char.zfill(str(fips[j]), width=5))\n            st_code = int((padstr)[0:2])\n            tfips = int(str(fips_to_code[st_code]) + padstr[2:])\n            cdf = tmin[(tmin['fips'] == tfips) & (tmin['year'] > 2000)]\n            cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n            cdf = cdf[(cdf != 0.00) & (cdf != -9.99)]\n            temp.append(np.mean(np.min(cdf, axis=1)))\n        ar = np.array(temp)\n        weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n        tmin_averaged.append(weighted_avg)\n\ntmin_std = []\nfor i in range(len(match_fips)):\n    if(not flag_multi_counties[i]):\n        fips = int(match_fips[i])\n        padstr = str(np.char.zfill(str(fips), width=5))\n        st_code = int((padstr)[0:2])\n        fips = int(str(fips_to_code[st_code]) + padstr[2:])\n        cdf = tmin[(tmin['fips'] == fips) & (tmin['year'] > 2000)]\n        cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float)  \n        cdf = cdf[(cdf != 0.00) & (cdf != -9.99)]\n        tmin_std.append(np.std(np.min(cdf, axis=1)))\n    else:\n        fips = match_fips[i]\n        temp = []\n        for j in range(len(fips)):\n            padstr = str(np.char.zfill(str(fips[j]), width=5))\n            st_code = int((padstr)[0:2])\n            tfips = int(str(fips_to_code[st_code]) + padstr[2:])\n            cdf = tmin[(tmin['fips'] == tfips) & (tmin['year'] > 2000)]\n            cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n            cdf = cdf[(cdf != 0.00) & (cdf != -9.99)]\n            temp.append(np.std(np.min(cdf, axis=1)))\n        ar = np.array(temp)\n        weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n        tmin_std.append(weighted_avg)\n        \ntmin_dt = []\nfor i in range(len(match_fips)):\n    if(not flag_multi_counties[i]):\n        fips = int(match_fips[i])\n        padstr = str(np.char.zfill(str(fips), width=5))\n        st_code = int((padstr)[0:2])\n        fips = int(str(fips_to_code[st_code]) + padstr[2:])\n        cdf = tmin[(tmin['fips'] == fips) & (tmin['year'] > 1980)]\n        cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float)  \n        cdf = cdf[(cdf != 0.00) & (cdf != -9.99)]\n        tmin_dt.append(np.mean(np.diff(np.max(cdf, axis=1))))\n    else:\n        fips = match_fips[i]\n        temp = []\n        for j in range(len(fips)):\n            padstr = str(np.char.zfill(str(fips[j]), width=5))\n            st_code = int((padstr)[0:2])\n            tfips = int(str(fips_to_code[st_code]) + padstr[2:])\n            cdf = tmin[(tmin['fips'] == tfips) & (tmin['year'] > 1980)]\n            cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n            cdf = cdf[(cdf != 0.00) & (cdf != -9.99)]\n            temp.append(np.mean(np.diff(np.max(cdf, axis=1))))\n        ar = np.array(temp)\n        weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n        tmin_dt.append(weighted_avg)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tmin_averaged_ranked = standardize_rank(tmin_averaged, 1)\ntmin_std_ranked = standardize_rank(tmin_std, 1)\ntmin_dt_ranked = standardize_rank(tmin_dt, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"occupation = pd.read_excel(path + \"Occupation_Data_OECS/MSA_M2019_dl.xlsx\", sheet_name=\"All May 2019 Data\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_locs_construction = []\nfor i in range(len(MSA)):\n    o = (occupation[occupation[\"area\"] == MSA[i]])\n    locs = o[o['occ_code'].str.startswith(\"47\")][\"loc_quotient\"]\n    locs = np.array([float(locs.iloc[x]) for x in range(len(locs)) if type(locs.iloc[x]) is float])\n    all_locs_construction.append(np.nanmean(locs))\n\nall_locs_transportation = []\nfor i in range(len(MSA)):\n    o = (occupation[occupation[\"area\"] == MSA[i]])\n    locs = o[o['occ_code'].str.startswith(\"53\")][\"loc_quotient\"]\n    locs = np.array([float(locs.iloc[x]) for x in range(len(locs)) if type(locs.iloc[x]) is float])\n    all_locs_transportation.append(np.nanmean(locs))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_locs_construction_ranked = standardize_rank(all_locs_construction, 1)\nall_locs_transportation_ranked = standardize_rank(all_locs_transportation, 1)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nri = pd.read_csv(path + \"NRI_Table_Counties.csv\")\n\ndef nri_select(field):\n    res = []\n    for i in range(len(match_fips)):\n        if(not flag_multi_counties[i]):\n            fips = int(match_fips[i])\n            cdf = nri[nri[\"STCOFIPS\"] == fips]\n            try:\n                res.append(float(cdf[field]))\n            except:\n                res.append(np.nan)\n        else:\n            fips = match_fips[i]\n            temp = []\n            for j in range(len(fips)):\n                cdf = nri[nri[\"STCOFIPS\"] == fips[j]]\n                try:\n                    temp.append(float(cdf[field]))\n                except:\n                    temp.append(np.nan)\n            ar = np.array(temp)\n            weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n            res.append(weighted_avg)\n    return np.array(res)\n\nwinter_weather_risk = nri_select(\"WNTW_EALS\")\nicestorm_risk = nri_select(\"ISTM_EALS\")\nhail_risk = nri_select(\"HAIL_EALS\")\ncold_wave_risk = nri_select(\"CWAV_EALS\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"winter_weather_risk_ranked = standardize_rank(winter_weather_risk, 1)\nicestorm_risk_ranked = standardize_rank(icestorm_risk, 1)\nhail_risk_ranked = standardize_rank(hail_risk, 1)\ncold_wave_risk_ranked = standardize_rank(cold_wave_risk, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vulnerability_raw = vulnerability.copy()\nvulnerability_ranked = vulnerability.copy()\nvulnerability_aggr = vulnerability.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### EXPOSURE\n\nvulnerability_raw[\"score_cold_days\"] = score_cold_days\nvulnerability_raw[\"score_cold_wave\"] = score_cold_wave\nvulnerability_raw[\"score_heavy_snow\"] = score_heavy_snow\nvulnerability_raw[\"score_winter_conditions\"] = score_winter_conditions\n\nvulnerability_raw[\"fema_cold_events\"] = fema_cold_events\nvulnerability_raw[\"tmin_averaged\"] = tmin_averaged\nvulnerability_raw[\"tmin_std\"] = tmin_std\nvulnerability_raw[\"tmin_dt\"] = tmin_dt\n\n\n### SENSITIVITY\n\nvulnerability_raw[\"LOC_construction\"] = all_locs_construction\nvulnerability_raw[\"LOC_transportation\"] = all_locs_transportation\nvulnerability_raw[\"winter_weather_expected_annual_loss\"] = winter_weather_risk\nvulnerability_raw[\"icestorm_expected_annual_loss\"] = icestorm_risk\nvulnerability_raw[\"hail_risk_expected_annual_loss\"] = hail_risk\nvulnerability_raw[\"cold_wave_expected_annual_loss\"] = cold_wave_risk\nvulnerability_raw.to_csv(\"vulnerability_cities_us_cold_kpis_raw.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### EXPOSURE\n\nvulnerability_ranked[\"cold_days_ranked\"] = cold_days_ranked\nvulnerability_ranked[\"cold_wave_ranked\"] = cold_wave_ranked\nvulnerability_ranked[\"heavy_snow_ranked\"] = heavy_snow_ranked\nvulnerability_ranked[\"winter_conditions_ranked\"] = winter_conditions_ranked\n\nvulnerability_ranked[\"fema_cold_events_ranked\"] = fema_cold_events_ranked\nvulnerability_ranked[\"tmin_averaged_ranked\"] = tmin_averaged_ranked\nvulnerability_ranked[\"tmin_std_ranked\"] = tmin_std_ranked\nvulnerability_ranked[\"tmin_dt_ranked\"] = tmin_dt_ranked\n\n### SENSITIVITY\n\nvulnerability_ranked[\"LOC_construction_ranked\"] = all_locs_construction_ranked\nvulnerability_ranked[\"LOC_transportation_ranked\"] = all_locs_transportation_ranked\nvulnerability_ranked[\"winter_weather_expected_annual_loss_ranked\"] = winter_weather_risk_ranked\nvulnerability_ranked[\"icestorm_expected_annual_loss_ranked\"] = icestorm_risk_ranked\nvulnerability_ranked[\"hail_risk_expected_annual_loss_ranked\"] = hail_risk_ranked\nvulnerability_ranked[\"cold_wave_expected_annual_loss_ranked\"] = cold_wave_risk_ranked\nvulnerability_ranked.to_csv(\"vulnerability_cities_us_cold_kpis_ranked.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def aggr_and_rescale(arr):\n    avg_arr = np.nanmean(arr, axis=1)\n    rescaled_avg_arr = standardize_rank(avg_arr)\n    return rescaled_avg_arr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_expo_cold_kpis = np.c_[cold_days_ranked, cold_wave_ranked, heavy_snow_ranked, winter_conditions_ranked,\n                          fema_cold_events_ranked,\n                          tmin_averaged_ranked, tmin_std_ranked, tmin_dt_ranked]\n\naggr_expo_cold_kpis = aggr_and_rescale(all_expo_cold_kpis)\n\nall_sens_cold_kpis = np.c_[all_locs_construction_ranked, all_locs_transportation_ranked, winter_weather_risk_ranked, \n                           icestorm_risk_ranked, hail_risk_ranked, cold_wave_risk_ranked]\n\naggr_sens_cold_kpis = aggr_and_rescale(all_sens_cold_kpis)\n\naggr_cold_kpi = aggr_and_rescale(np.c_[aggr_expo_cold_kpis, aggr_sens_cold_kpis])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vulnerability_aggr = vulnerability.copy()\n\nvulnerability_aggr[\"aggr_cold_kpis\"] = aggr_cold_kpi\nvulnerability_aggr[\"aggr_exposure_cold_kpis\"] = aggr_expo_cold_kpis\nvulnerability_aggr[\"aggr_sensitivity_cold_kpis\"] = aggr_sens_cold_kpis\n\nvulnerability_aggr.to_csv(\"vulnerability_cities_us_cold_aggr_kpis.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}