{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.core.display import HTML\nHTML(\"\"\"\n\n\n\"\"\")","execution_count":null,"outputs":[]},{"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":{},"cell_type":"markdown","source":"| severe_storm_score | Calculated score based on survey response | CDP Survey |\n|---------------------------|---------------------------------------------------------------------------|--------------------------|\n| fema_extreme_events | # of FEMA severe storm disasters | FEMA National Risk Index |\n| percent_built_before_1999 | Percentage of buildings built before 1999 | US Census |\n| precip_averaged | Minimum temp per month averaged over the last 20 years | NOAA tmin data |\n| precip_std | Variability in minimum temp per month averaged over the last 20 years | NOAA tmin data |\n| precip_dt | Rate of change in minimum temp per month averaged over the last 20 years | NOAA tmin data |"},{"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":"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","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":"scores_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}\n\ndef 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"risk_severe_wind, score_severe_wind = return_risk_and_score_per_county(\"Storm and wind > Severe wind\")\n\nrisk_tornado, score_tornado = return_risk_and_score_per_county(\"Storm and wind > Tornado\")\n\nrisk_cyclone, score_cyclone = return_risk_and_score_per_county(\"Storm and wind > Cyclone (Hurricane / Typhoon)\")\n\nrisk_tropical, score_tropical = return_risk_and_score_per_county(\"Storm and wind > Tropical storm\")\n\nrisk_thunderstorm, score_thunderstorm = return_risk_and_score_per_county(\"Storm and wind > Lightning / thunderstorm\")\n\nrisk_storm_surge, score_storm_surge = return_risk_and_score_per_county(\"Storm and wind > Storm surge\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"severe_storm_score = score_severe_wind + score_tornado + score_cyclone + score_tropical + score_thunderstorm + score_storm_surge\nsevere_storm_ranked = standardize_rank(severe_storm_score, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fema = pd.read_csv(\"../input/femadisasters/DisasterDeclarationsSummaries.csv\")\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_extreme_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\"] == \"Typhoon\") | (fema[\"incidentType\"] == \"Severe Storm(s)\") | (fema[\"incidentType\"] == \"Tornado\")| (fema[\"incidentType\"] == \"Hurricane\") ]\n cdf = cdf[cdf[\"fips\"] == fips]\n fema_extreme_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\"] == \"Typhoon\") | (fema[\"incidentType\"] == \"Severe Storm(s)\") | (fema[\"incidentType\"] == \"Tornado\")| (fema[\"incidentType\"] == \"Hurricane\") ] \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_extreme_events.append(weighted_avg)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fema_extreme_events_ranked = standardize_rank(fema_extreme_events, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building = pd.read_csv(\"../input/housingyear/Year_Housing/Building_Year.csv\", skiprows=1)\nbuilding.head()\nbuilding = building.loc[:, ~building.columns.str.startswith('Margin')]\nbuilding[\"id\"] = building[\"id\"].str[9:].astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"percent_built_before_1999 = []\nfor i in range(len(match_fips)):\n if(not flag_multi_counties[i]):\n fips = int(match_fips[i])\n cdf = building[building[\"id\"] == fips]\n percent_built_before_1999.append(np.sum(np.array(cdf)[0][6:])/(np.array(cdf)[0][2]))\n else:\n fips = match_fips[i]\n temp = []\n for j in range(len(fips)):\n cdf = building[building[\"id\"] == fips[j]]\n temp.append(np.sum(np.array(cdf)[0][6:])/(np.array(cdf)[0][2]))\n ar = np.array(temp)\n weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n percent_built_before_1999.append(weighted_avg)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"percent_built_before_1999_ranked = standardize_rank(percent_built_before_1999, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"str_precip = open('../input/precipitationyearnoaa/climdiv-pcpncy-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_precip), header=None, chunksize=500000, names=columns,\n converters={h:str for h in columns}):\n precip = df\n\nyears = []\nfips = []\nfor i in range(len(precip)):\n fips.append(int(str(precip[\"Code\"].iloc[i])[0:5]))\n years.append(int(str(precip[\"Code\"].iloc[i])[7:11]))\n \nprecip[\"year\"] = np.array(years).astype(int)\nprecip[\"fips\"] = np.array(fips).astype(int)\n\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\n\nprecip_max_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 = precip[(precip['fips'] == fips) & (precip['year'] > 2000)]\n cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float)\n cdf = cdf[(cdf != -9.99)]\n precip_max_averaged.append(np.mean(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 = precip[(precip['fips'] == tfips) & (precip['year'] > 2000)]\n cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n cdf = cdf[(cdf != -9.99)]\n temp.append(np.mean(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 precip_max_averaged.append(weighted_avg)\n\nprecip_max_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 = precip[(precip['fips'] == fips) & (precip['year'] > 2000)]\n cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n cdf = cdf[(cdf != -9.99)]\n precip_max_std.append(np.std(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 = precip[(precip['fips'] == tfips) & (precip['year'] > 2000)]\n cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n cdf = cdf[(cdf != -9.99)]\n temp.append(np.mean(np.std(cdf, axis=1)))\n ar = np.array(temp)\n weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n precip_max_std.append(weighted_avg)\n\nprecip_max_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 = precip[(precip['fips'] == fips) & (precip['year'] > 1980)]\n cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n cdf = cdf[(cdf != -9.99)]\n precip_max_dt.append(np.median(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 = precip[(precip['fips'] == tfips) & (precip['year'] > 1980)]\n cdf = cdf.drop(columns=[\"Code\", \"year\", \"fips\"]).astype(np.float) \n cdf = cdf[(cdf != -9.99)]\n temp.append(np.median(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 precip_max_dt.append(weighted_avg)\n\nprecip_max_averaged_ranked = standardize_rank(precip_max_averaged, 1)\nprecip_max_std_ranked = standardize_rank(precip_max_std, 1)\nprecip_max_dt_ranked = standardize_rank(precip_max_dt, 1)","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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"avalanche_eals = nri_select(\"AVLN_EALS\")\nearthquake_eals = nri_select(\"ERQK_EALS\")\nstrong_wind_eals = nri_select(\"SWND_EALS\")\ntornado_eals = nri_select(\"TRND_EALS\")\n\navalanche_eals_ranked = standardize_rank(avalanche_eals, 1)\nearthquake_eals_ranked = standardize_rank(earthquake_eals, 1)\nstrong_wind_eals_ranked = standardize_rank(strong_wind_eals, 1)\ntornado_eals_ranked = standardize_rank(tornado_eals, 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_severe_storm\"] = severe_storm_score\nvulnerability_raw[\"fema_extreme_events\"] = fema_extreme_events\nvulnerability_raw[\"precip_max_averaged\"] = precip_max_averaged\nvulnerability_raw[\"precip_max_std\"] = precip_max_std\nvulnerability_raw[\"precip_max_dt\"] = precip_max_dt\n\n### ECONOMIC\nvulnerability_raw[\"percent_built_before_1999\"] = percent_built_before_1999\nvulnerability_raw[\"avalanche_expected_annual_loss\"] = avalanche_eals\n\nvulnerability_raw[\"earthquake_expected_annual_loss\"] = earthquake_eals\nvulnerability_raw[\"strong_wind_expected_annual_loss\"] = strong_wind_eals\nvulnerability_raw[\"tornado_expected_annual_loss\"] = tornado_eals\n\nvulnerability_raw.to_csv(\"vulnerability_cities_us_severe_storms_kpis_raw.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### EXPOSURE \n\nvulnerability_ranked[\"severe_storm_ranked\"] = severe_storm_ranked\nvulnerability_ranked[\"fema_extreme_events_ranked\"] = fema_extreme_events_ranked\nvulnerability_ranked[\"precip_max_averaged_ranked\"] = precip_max_averaged_ranked\nvulnerability_ranked[\"precip_max_std_ranked\"] = precip_max_std_ranked\nvulnerability_ranked[\"precip_max_dt_ranked\"] = precip_max_dt_ranked\n\n\n### ECONOMIC\nvulnerability_ranked[\"percent_built_before_1999_ranked\"] = percent_built_before_1999_ranked\nvulnerability_ranked[\"avalanche_expected_annual_loss_ranked\"] = avalanche_eals_ranked\n\nvulnerability_ranked[\"earthquake_expected_annual_loss_ranked\"] = earthquake_eals_ranked\nvulnerability_ranked[\"strong_wind_expected_annual_loss_ranked\"] = strong_wind_eals_ranked\nvulnerability_ranked[\"tornado_expected_annual_loss_ranked\"] = tornado_eals_ranked\n\nvulnerability_ranked.to_csv(\"vulnerability_cities_us_severe_storms_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_severe_kpis = np.c_[severe_storm_ranked, fema_extreme_events_ranked, precip_max_averaged_ranked, precip_max_std_ranked,\n precip_max_dt_ranked]\n\naggr_expo_severe_kpis = aggr_and_rescale(all_expo_severe_kpis)\n\nall_sens_severe_kpis = np.c_[percent_built_before_1999_ranked, avalanche_eals_ranked, earthquake_eals_ranked, strong_wind_eals_ranked, \n tornado_eals_ranked]\n\naggr_sens_severe_kpis = aggr_and_rescale(all_sens_severe_kpis)\n\naggr_severe_kpi = aggr_and_rescale(np.c_[aggr_expo_severe_kpis, aggr_sens_severe_kpis])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vulnerability_aggr = vulnerability.copy()\n\nvulnerability_aggr[\"aggr_severe_kpis\"] = aggr_severe_kpi\nvulnerability_aggr[\"aggr_exposure_severe_kpis\"] = aggr_expo_severe_kpis\nvulnerability_aggr[\"aggr_sensitivity_severe_kpis\"] = aggr_sens_severe_kpis\n\nvulnerability_aggr.to_csv(\"vulnerability_cities_us_severe_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}