{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.core.display import HTML\nHTML(\"\"\"\n\n\n\"\"\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Deriving all vulnerability KPIs for cities in the US\n\n## This notebook uses data generated from the different vulnerability_* notebooks. \n\nDescription of KPIs: (the climate hazard-specific ones can be found in the different notebooks)\n\n| KPI Name | KPI Interpretation | Public Data Source |\n|----------------------|---------------------------------------------------------------|--------------------------|\n| human_pressure | Human pressure on environment footprint | SEDAC |\n| aqi_mean | Mean of air quality index over last three years | Kaggle Dataset |\n| daqi | Change in air quality index year-by-year | Kaggle Dataset |\n| air_pollution_trends | Trend in air pollution (PM2.5) | County Health Rankings |\n| pop_density | population density | Survey |\n| emissions_per_capita | emissions per capita in different counties | Jones & Kammen (2014) |\n| emissions_per_km2 | emissions per land area in different counties | Jones & Kammen (2014) |\n| expected_annual_loss | overall expected annual loss from different climate disasters | FEMA National Risk Index |\n| risk_city_sentiment | city sentiment on the risks they face | CDP Survey |"},{"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\n\nfrom tqdm import *\nfrom pyproj import Transformer\nimport ast\nfrom sklearn.preprocessing import RobustScaler, QuantileTransformer, StandardScaler\n","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/\" ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vulnerability = 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":{},"cell_type":"markdown","source":"## First, we use SEDAC data to calculate the human footprint on the environment"},{"metadata":{"trusted":true},"cell_type":"code","source":"fp = path + 'Human_Footprint/human_footprint.tif'\nimg = rasterio.open(fp)\nimg_array = img.read(1)\ntransformer = Transformer.from_crs(\"epsg:4326\", img.crs)\n\ndef get_value(raster, array, lon, lat):\n coords = transformer.transform(lat, lon) \n py, px = raster.index(coords[0], coords[1])\n return array[py, px]\n\nplt.imshow(img_array, norm=colors.LogNorm())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"human_pressure = []\nfor i in tqdm(range(len(vulnerability))):\n coords = ((vulnerability['loc'].iloc[i].replace(\"POINT (\", \"\")).replace(\")\",\"\")).split()\n lon = float(coords[0])\n lat = float(coords[1])\n human_pressure.append(get_value(img, img_array, lon, lat))\n \ndirection_human_pressure = 1.0\nhuman_pressure_score = standardize_rank(np.array(human_pressure), direction=direction_human_pressure)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's calculate some emission / air pollution statistics. "},{"metadata":{"trusted":true},"cell_type":"code","source":"aqi = pd.read_csv(path + 'AQI_by_year.csv')\naccts_2020 = list(vulnerability[\"Account Number\"])\n\naqi_mean = []\nfor i in range(len(accts_2020)):\n try:\n aa = (aqi[aqi[\"account_number\"] == accts_2020[i]]).sort_values(\"year\")[\"arithmetic_mean\"]\n da = np.mean(aa.iloc[-3:])\n aqi_mean.append(da)\n except:\n aqi_mean.append(np.nan)\n \ndaqi = []\nfor i in range(len(accts_2020)):\n try:\n aa = (aqi[aqi[\"account_number\"] == accts_2020[i]]).sort_values(\"year\")[\"arithmetic_mean\"]\n da = np.mean(np.diff(aa)/aa.iloc[0])\n daqi.append(da)\n except:\n daqi.append(np.nan)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"direction_air_quality = 1.0\naqi_score = standardize_rank(np.array(aqi_mean), direction=direction_air_quality)\n\ndirection_air_quality_trend = 1.0\ndaqi_score = standardize_rank(np.array(daqi), direction=direction_air_quality_trend)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"trends = pd.read_csv(path + \"CHR_trends_csv_2020.csv\", encoding='latin-1')\nct = np.char.zfill(np.array(trends['countycode']).astype(str), 3)\nst = np.char.zfill(np.array(trends['statecode']).astype(str), 2)\ntrends['fips'] = np.core.defchararray.add(st, ct).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"air_pollution_trends = []\nfor i in range(len(match_fips)):\n if(not flag_multi_counties[i]):\n fips = match_fips[i]\n cdf = trends[(trends['measurename'] == \"Air pollution - particulate matter\") & (trends['fips'] == fips)]\n rv = np.array(cdf['rawvalue']).astype(float)\n rvd = (np.mean(rv[-3:]) - np.mean(rv[0:3]))/(np.mean(rv[0:3]))\n air_pollution_trends.append(np.float(rvd))\n else:\n fips = match_fips[i]\n temp = []\n for j in range(len(fips)):\n cdf = trends[(trends['measurename'] == \"Air pollution - particulate matter\") & (trends['fips'] == fips[j])]\n rv = np.array(cdf['rawvalue']).astype(float)\n rvd = (np.mean(rv[-3:]) - np.mean(rv[0:3]))/(np.mean(rv[0:3]))\n temp.append(np.float(rvd))\n ar = np.array(temp)\n weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n air_pollution_trends.append(weighted_avg)\n\ndirection_air_pollution_trend = 1.0\nair_pollution_trend_score = standardize_rank(np.array(air_pollution_trends), direction=direction_air_pollution_trend)","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":"population = np.array([float(fetch_answer_individual(accts_2020[i], '0.5', column_number=1.0)[\"Response Answer\"]) for i in range(len(accts_2020))])\nareas_km2 = np.array([float(fetch_answer_individual(accts_2020[i], '0.6')[\"Response Answer\"]) for i in range(len(accts_2020))])\nareas_km2[cities == \"Greenbelt, MD\"] /= 10**3 # Greenbelt reported in the wrong units\npop_density = population/(areas_km2)\ndirection_pop_density = 1\nlog_pop_density_score = standardize_rank(np.log10(pop_density), direction=direction_pop_density)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"carbon_footprint = pd.read_excel(\"/kaggle/input/joneskammen/Jones-Kammen-2014-Zip-City-County-Results.xlsx\", sheet_name=\"County results\")\ncarbon_footprint[\"County\"] = carbon_footprint[\"County\"].str.title()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emissions_per_capita = []\nfor i in range(len(counties)):\n if(not flag_multi_counties[i]):\n ct = ((counties[i].replace(\" County\", \"\")).replace(\"'\", \"\").replace(\".\",\"\").replace(\" Parish\", \"\")).title()\n cdf = carbon_footprint[(carbon_footprint[\"County\"].str[:-1] == ct) & (carbon_footprint[\"State\"] == states[i])]\n emissions_per_capita.append(float(cdf[\"Total County Carbon Footprint (tCO2e/yr)\"]/cdf[\"Population\"]))\n else:\n fips = ast.literal_eval(counties[i])\n temp = []\n for j in range(len(fips)):\n ct = ((fips[j].replace(\" County\", \"\")).replace(\"'\", \"\").replace(\".\",\"\").replace(\" Parish\", \"\")).title()\n cdf = carbon_footprint[(carbon_footprint[\"County\"].str[:-1] == ct) & (carbon_footprint[\"State\"] == states[i])]\n temp.append((float(cdf[\"Total County Carbon Footprint (tCO2e/yr)\"]/cdf[\"Population\"])))\n ar = np.array(temp)\n weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n emissions_per_capita.append(weighted_avg)\n ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emissions_per_km2 = []\nfor i in range(len(counties)):\n if(not flag_multi_counties[i]):\n ct = ((counties[i].replace(\" County\", \"\")).replace(\"'\", \"\").replace(\".\",\"\").replace(\" Parish\", \"\")).title()\n cdf = carbon_footprint[(carbon_footprint[\"County\"].str[:-1] == ct) & (carbon_footprint[\"State\"] == states[i])]\n emissions_per_km2.append(float(cdf[\"Total County Carbon Footprint (tCO2e/yr)\"]/areas_km2[i]))\n else:\n fips = ast.literal_eval(counties[i])\n temp = []\n for j in range(len(fips)):\n ct = ((fips[j].replace(\" County\", \"\")).replace(\"'\", \"\").replace(\".\",\"\").replace(\" Parish\", \"\")).title()\n cdf = carbon_footprint[(carbon_footprint[\"County\"].str[:-1] == ct) & (carbon_footprint[\"State\"] == states[i])]\n temp.append((float(cdf[\"Total County Carbon Footprint (tCO2e/yr)\"]/areas_km2[i])))\n ar = np.array(temp)\n weighted_avg = np.average(ar, weights=np.array(ast.literal_eval(counties_weights[i])))\n emissions_per_km2.append(weighted_avg)\n ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"direction_emissions_per_capita = 1\nemissions_per_capita_score = standardize_rank(np.array(emissions_per_capita), direction=direction_emissions_per_capita)\n\ndirection_emissions_per_km2 = 1\nemissions_per_km2_score = standardize_rank(np.log10(emissions_per_km2), direction=direction_emissions_per_km2)","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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"expected_annual_loss = nri_select(\"EAL_SCORE\")\ndirection_expected_annual_loss = 1\nexpected_annual_loss_score = standardize_rank(expected_annual_loss, direction=direction_expected_annual_loss)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from nltk.sentiment.vader import SentimentIntensityAnalyzer\nfrom nltk import tokenize\n\nsid = SentimentIntensityAnalyzer()\n\ndef polarity_score_city(question, column, sent_thresh=2):\n subset = cities_2020[cities_2020['Question Number'] == question]\n subset_city_answers = subset[subset[\"Column Number\"] == column]\n \n compound_score = []\n for i in range(len(acct_number)):\n ss = list(subset_city_answers[subset_city_answers[\"Account Number\"] == acct_number[i]][\"Response Answer\"])\n ss = ' '.join(map(str, ss))\n sents = tokenize.sent_tokenize(ss)\n if(len(sents) > sent_thresh):\n cc = 0\n for s in sents:\n r = sid.polarity_scores(s)\n cc += r['compound']\n avg_cc = cc/len(sents)\n compound_score.append(avg_cc)\n else:\n compound_score.append(np.nan)\n return (compound_score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"risk_city_sentiment = polarity_score_city(\"2.1\", 12.0)\nrisk_city_sentiment_ranked = standardize_rank(risk_city_sentiment, -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[\"log_pop_density\"] = np.log10(pop_density)\nvulnerability_raw[\"human_pressure\"] = human_pressure\nvulnerability_raw[\"expected_annual_loss\"] = expected_annual_loss\nvulnerability_raw[\"risk_city_sentiment\"] = risk_city_sentiment\n\n### EMISSIONS & AIR POLLUTION\n\nvulnerability_raw[\"aqi_mean\"] = aqi_mean\nvulnerability_raw[\"daqi\"] = daqi\nvulnerability_raw[\"air_pollution_trend\"] = air_pollution_trends\nvulnerability_raw[\"emissions_per_capita\"] = emissions_per_capita\nvulnerability_raw[\"emissions_per_km2\"] = emissions_per_km2\n\nvulnerability_raw.to_csv(\"vulnerability_cities_us_general_kpis_raw.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### EXPOSURE\n\nvulnerability_ranked[\"log_pop_density_score\"] = log_pop_density_score\nvulnerability_ranked[\"human_pressure_score\"] = human_pressure_score\nvulnerability_ranked[\"expected_annual_loss_score\"] = expected_annual_loss_score\nvulnerability_ranked[\"risk_city_sentiment_score\"] = risk_city_sentiment_ranked\n\n### EMISSIONS & AIR POLLUTION\n\nvulnerability_ranked[\"aqi_score\"] = aqi_score\nvulnerability_ranked[\"daqi_score\"] = daqi_score\nvulnerability_ranked[\"air_pollution_trend_score\"] = air_pollution_trend_score\nvulnerability_ranked[\"emissions_per_capita_score\"] = emissions_per_capita_score\nvulnerability_ranked[\"emissions_per_km2_score\"] = emissions_per_km2_score\nvulnerability_ranked.to_csv(\"vulnerability_cities_us_general_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 = (avg_arr - np.nanmin(avg_arr))/(np.nanmax(avg_arr) - np.nanmin(avg_arr))\n return rescaled_avg_arr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_vul_general_kpis = np.c_[log_pop_density_score, human_pressure_score, expected_annual_loss_score, \n risk_city_sentiment_ranked, aqi_score, daqi_score, air_pollution_trend_score, emissions_per_capita_score,\n emissions_per_km2_score]\n\naggr_vul_general_kpis = aggr_and_rescale(all_vul_general_kpis)\n\nvulnerability_aggr[\"aggr_vul_general_kpis\"] = aggr_vul_general_kpis\nvulnerability_aggr.to_csv(\"vulnerability_cities_us_general_aggr_kpis.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vul_aggr_heat = pd.read_csv(\"../input/derivedkpis/KPIs/vulnerability_cities_us_heat_aggr_kpis.csv\")\nvul_aggr_cold = pd.read_csv(\"../input/derivedkpis/KPIs/vulnerability_cities_us_cold_aggr_kpis.csv\")\nvul_aggr_drought = pd.read_csv(\"../input/derivedkpis/KPIs/vulnerability_cities_us_drought_aggr_kpis.csv\")\nvul_aggr_flood = pd.read_csv(\"../input/derivedkpis/KPIs/vulnerability_cities_us_floods_aggr_kpis.csv\")\nvul_aggr_severe = pd.read_csv(\"../input/derivedkpis/KPIs/vulnerability_cities_us_severe_aggr_kpis.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggr_vul_general_kpis = aggr_and_rescale(all_vul_general_kpis)\naggr_vul_all = np.c_[aggr_vul_general_kpis, vul_aggr_heat[\"aggr_heat_kpis\"], vul_aggr_cold[\"aggr_cold_kpis\"], vul_aggr_drought[\"aggr_drought_kpis\"]\n , vul_aggr_flood[\"aggr_floods_kpis\"], vul_aggr_severe[\"aggr_severe_kpis\"]]\n\naggr_vul_all_kpis = aggr_and_rescale(aggr_vul_all)\n\nvulnerability_aggr_all = vulnerability.copy()\n\nvulnerability_aggr_all[\"aggr_heat_kpis\"] = vul_aggr_heat[\"aggr_heat_kpis\"]\nvulnerability_aggr_all[\"aggr_cold_kpis\"] = vul_aggr_cold[\"aggr_cold_kpis\"]\nvulnerability_aggr_all[\"aggr_drought_kpis\"] = vul_aggr_drought[\"aggr_drought_kpis\"]\nvulnerability_aggr_all[\"aggr_flood_kpis\"] = vul_aggr_flood[\"aggr_floods_kpis\"]\nvulnerability_aggr_all[\"aggr_severe_kpis\"] = vul_aggr_severe[\"aggr_severe_kpis\"]\nvulnerability_aggr_all[\"aggr_vul_all_kpis\"] = aggr_vul_all_kpis\n\nvulnerability_aggr_all.to_csv(\"vulnerability_cities_us_all_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}