library(sf)

library(data.table)

countries = st_read("https://michaelminn.net/tutorials/data/2024-natural-earth-countries.geojson")

countries = countries[,c("NAME", "ISO_A3", "geometry")]

names(countries) = c("Country_Name", "Country_Code", "geometry")

indicators = c(
"Population" = 			"SP.POP.TOTL",
"Working_Age_Percent" =		"SP.POP.1564.TO.ZS",
"Urban_Population_Percent" = 	"SP.URB.TOTL.IN.ZS",
"Rural_Population_Percent" = 	"SP.RUR.TOTL.ZS",
"Life_Expectancy_Years" = 	"SP.DYN.LE00.IN",
"Percent_Secondary_School" = 	"SE.SEC.CUAT.UP.ZS",
"Percent_Bachelors" =		"SE.TER.CUAT.BA.ZS",
"Land_Area_Sq_km" = 		"AG.LND.TOTL.K2",
"Forest_Area_Percent" = 	"AG.LND.FRST.ZS",
"Agricultural_Land_Percent" = 	"AG.LND.AGRI.ZS",
"Arable_Land_Percent" = 	"AG.LND.ARBL.ZS",
"Low_Elevation_Percent" = 	"EN.POP.EL5M.ZS",
"Water_Percent" = 		"ER.H2O.FWTL.ZS",
"Irrigation_Percent" = 		"AG.LND.IRIG.AG.ZS",
"Fertilizer_Intensity" = 	"AG.CON.FERT.ZS",
"GDP_per_Capita_PPP" = 		"NY.GDP.PCAP.PP.CD",
"Gini_Index" = 			"SI.POV.GINI",
"Agriculture_Percent_GDP" = 	"NV.AGR.TOTL.ZS",
"Industry_Percent_GDP" = 	"NV.IND.TOTL.ZS",
"Military_Percent_GDP" = 	"MS.MIL.XPND.GD.ZS",
"Resource_Percent_GDP" =	"NY.GDP.TOTL.RT.ZS",
"Energy_MM_BTU_per_Capita" = 	"EG.USE.PCAP.KG.OE",
"CO2_Tonnes_per_Capita" = 	"EN.GHG.CO2.PC.CE.AR5",
"Renewable_Percent" = 		"EG.FEC.RNEW.ZS",
"Internet_Percent" = 		"IT.NET.USER.ZS",
"Male_Youth_Unemployment" = 	"SL.UEM.1524.MA.ZS",
"Attended_Births_Percent" = 	"SH.STA.BRTC.ZS",
"Contraceptive_Percent" = 	"SP.DYN.CONU.ZS",
"Diabetes_Percent" = 		"SH.STA.DIAB.ZS",
"HIV_Prevalence_Percent" = 	"SH.DYN.AIDS.ZS",
"Immunized_DPT_Percent" = 	"SH.IMM.IDPT",
"Immunized_Measles_Percent" = 	"SH.IMM.MEAS",
"Mortality_Communicable_Percent" = "SH.DTH.COMM.ZS",
"Mortality_Infant_per_1K" = 	"SP.DYN.IMRT.IN",
"Mortality_Injury_Percent" = 	"SH.DTH.INJR.ZS",
"Mortality_Maternal_per_100K" = "SH.STA.MMRT",
"Mortality_Non-Communicable_Deaths_Percent" = "SH.DTH.NCOM.ZS",
"Mortality_Traffic_per_100K" = 	"SH.STA.TRAF.P5",
"Surgeons_per_100K" = 		"SH.MED.SAOP.P5",
"Teen_Fertility_per_1K" = 	"SP.ADO.TFRT",
"Tuberculosis_per_100K" = 	"SH.TBS.INCD",
"Undernourished_Percent" = 	"SN.ITK.DEFC.ZS")


for (varname in names(indicators)) {
	indicator = indicators[varname]

	print(paste(varname, indicator))

	url = paste0("https://api.worldbank.org/v2/en/indicator/", indicator, "?downloadformat=csv")

	download.file(url, "temp.zip")

	name_list = unzip("temp.zip", list=T)$Name

	name_list = name_list[!grepl("Metadata", name_list)]

	data = read.csv(unz("temp.zip", name_list), skip=4)

	data[,varname] = sapply(1:nrow(data), function(x) suppressWarnings(as.numeric(data[x, last(which(!is.na(data[x,])))])))

	data = data[,c("Country.Code", varname)]

	names(data) = c("Country_Code", varname)

	countries = merge(countries, data)
}

countries = countries[order(countries$Country_Name),]

# kg oil to MM BTU

countries[, "Energy_MM_BTU_per_Capita"] = countries[, "Energy_MM_BTU_per_Capita"] * 39653 / 1000000

summary(countries)

st_write(countries, "2024-world-bank-indicators.geojson", delete_dsn=T)

data = st_drop_geometry(countries)

write.csv(data, "2024-world-bank-indicators.csv", na="", row.names=F)

projected = st_transform(countries, st_crs("ESRI:54030")) # World Robinson

plot(projected["Energy_MM_BTU_per_Capita"], logz=T)