![\"Cumulative](\"https://grantwitness.org/cms_media/noaa-obligations-curve-2026-09-25.png\"/)
Unlike some other agencies we track, NOAA grants lack an agency-specific data source for grants.\nFor NOAA, we rely primarily on USAspending, DOGE, and crowd-sourced information to make determinations on grant status.
\nDOGE provides data on terminated awards, although we have found that it is incomplete or occasionally incorrect, so we generally do manual review of DOGE data or use it only in combination with other sources.
\nWe use this source primarily for metadata including grantee information, grant assistance listing (CFDA), grant notice of funding announcement (NOFO), and total obligated and total outlaid amounts.\nWe also track changes to project end dates (i.e. project end date moved earlier, a \"cutoff\", or later, a \"reextension\"), and deobligation events (i.e. a negative amount obligated) which we use in combination as signals of disruption.
\nFile C is released roughly once a month (except that Sep/Oct are combined into one release).\nWe use this source primarily for outlay actions which we use to confirm reinstatements.
\nAn award is marked as terminated if it has been reported as such directly to us by a PI or via a trusted source, if it was listed by DOGE as being terminated, or if we see both a project end date cutoff and a deobligation in USAspending data.
\nAn award is labeled as reinstated if it has been reported as such directly to us by a PI, or if there is a re-extension of project end date or a positive obligation after a deobligation in USAspending awards data.\nReinstatements are confirmed by verifying that an award has received outlays since the reinstatement.
\nThe funds obligated and spent at time of earliest disruption are calculated using the current total obligated amount minus obligations since the earliest termination date and the total outlays minus outlays since the earliest termination date, respectively.\nThe outlay and obligation events come from File C, while the current total obligated and outlaid numbers are from the more regularly updating USAspending awards data.\nAs a result, these estimates may fluctuate week-to-week.
","doi":"https://doi.org/10.59350/74brw-em383","guid":"https://doi.org/10.59350/74brw-em383","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1790553600,"rid":"c2qz2-pwj80","summary":"Unlike some other agencies we track, NOAA grants lack an agency-specific data source for grants. For NOAA, we rely primarily on USAspending, DOGE, and crowd-sourced information to make determinations on grant status.","tags":["Methods","National Oceanic And Atmospheric Administration"],"title":"NOAA Grant Disruption Methodology","updated_at":1790621962,"url":"https://grantwitness.org/noaa/updates/methodology","version":"v1"}},{"document":{"authors":[{"contributor_roles":[],"family":"Shan","given":"Mridula","url":"https://orcid.org/0000-0003-4613-2669"},{"affiliation":[{"id":"https://ror.org/019jywm96","name":"rOpenSci"}],"contributor_roles":[],"family":"Ross","given":"Noam","url":"https://orcid.org/0000-0002-2136-0000"}],"blog":{"authors":null,"community_id":"3446712e-ab03-4bb7-8641-571dcda8d8cd","created":1787702400,"current_feed_url":null,"description":"Tracking disruptions to U.S. federal science funding, grant by grant.","doi":"https://doi.org/10.59350/grantwitness","favicon":"https://rogue-scholar.org/api/communities/3446712e-ab03-4bb7-8641-571dcda8d8cd/logo","feed_format":"application/atom+xml","feed_url":"https://grantwitness.org/updates/rogue-scholar.xml","filter":null,"generator":"Other","home_page_url":"https://grantwitness.org/updates?type=updates","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"grantwitness","status":"active","subfield":"3321","title":"Grant Witness","updated":1790553600,"use_api":null},"blog_name":"Grant Witness","blog_slug":"grantwitness","content_html":"At Grant Witness, we have been tracking the damage of funding terminations across federal agencies since March 2025. We have focused on our core research agencies (NIH and NSF) and expanded to other agencies in science and public health such as EPA and CDC. Yet scientific research, development, and application are core functions across the federal government, and we are always striving to expand our view to better understand how the current administration is upending the American science enterprise.
\nToday we're expanding both our termination and grantmaking trackers to include the National Oceanic and Atmospheric Administration (NOAA). NOAA is a scientific agency within the United States Department of Commerce. Its stated mission is \"to understand and predict our changing environment, from the deep sea to outer space, and to manage and conserve America's coastal marine resources.\" It includes the National Weather Service.
\nNOAA's work has been targeted from the start of the administration. The Trump administration has engaged in staff firings, suppression of scientific reports, and proposed budget cuts of more than 25% for both FY 2026 and 2027.
\nLike many other agencies, NOAA terminated research and support grants en masse starting in April 2025 under the disastrous tenure of the \"Department of Government Efficiency\" (DOGE). Neither NOAA nor DOGE has published a reliable accounting of these terminations. Grant Witness, through reporting and data mining, has identified 26 NOAA grants that have been terminated, only 4 of which have been reinstated.
\nThese grants span 13 government programs (CFDAs, Table\u00a01) across 4 of the 6 line offices at NOAA (Table\u00a02). Together, they represent a loss of $28,237,466 in funding to states, cities, nonprofits and universities. Climate research is the most affected area, but a variety of programs in coastal and marine management have been shuttered.
\n| CFDA Program Title | \nDisrupted Grants | \nAmt Remaining at Disruption | \n
|---|---|---|
| Climate and Atmospheric Research | \n11 | \n$11,217,919 | \n
| Cooperative Institute Inter Agency Funded Activities | \n2 | \n$1,927,901 | \n
| Marine Sanctuary Program | \n2 | \n$235,000 | \n
| Office for Coastal Management | \n2 | \n$1,470,012 | \n
| Chesapeake Bay Studies | \n1 | \n$40,700 | \n
| Coastal Zone Management Administration Awards | \n1 | \n$179,600 | \n
| Coastal Zone Management Estuarine Research Reserves | \n1 | \n$8,486.71 | \n
| Cooperative Science and Education Program | \n1 | \n$11,600 | \n
| Habitat Conservation | \n1 | \n$15.51 | \n
| National Oceanic and Atmospheric Administration Noaa Cooperative Institutes | \n1 | \n$2,757,930 | \n
| Noaa Mission Related Education Awards | \n1 | \n$364,360 | \n
| Pacific Fisheries Data Program | \n1 | \n$766,710 | \n
| Sea Grant Support | \n1 | \n$9,257,231 | \n
| Line Office Name | \nDisrupted Grants | \nAmt Remaining at Disruption | \n
|---|---|---|
| Office of Oceanic and Atmospheric Research | \n15 | \n$25,160,981 | \n
| National Ocean Service | \n6 | \n$1,893,099 | \n
| National Marine Fisheries Service | \n4 | \n$819,026 | \n
| Office of the Under Secretary | \n1 | \n$364,360 | \n
As one would expect, coastal states have borne the brunt of these terminations. New York has the greatest number of disrupted NOAA grants followed by Washington, Massachusetts, and New Jersey (tied at three disrupted grants) (Table\u00a03). The states with the greatest loss in funding include Washington, North Crolina, and New Jersey.
\n| State | \nDisrupted Grants | \nRemaining at Disruption | \n
|---|---|---|
| New York | \n4 | \n$787,780 | \n
| Washington | \n3 | \n$11,918,505 | \n
| New Jersey | \n3 | \n$4,685,831 | \n
| Massachusetts | \n3 | \n$384,447 | \n
| North Carolina | \n2 | \n$5,725,610 | \n
| Virginia | \n2 | \n$124,667 | \n
| California | \n1 | \n$1,921,064 | \n
| Louisiana | \n1 | \n$1,126,012 | \n
| Oregon | \n1 | \n$766,710 | \n
| Rhode Island | \n1 | \n$344,000 | \n
| Wyoming | \n1 | \n$271,297 | \n
| Michigan | \n1 | \n$89,519 | \n
| Illinois | \n1 | \n$75,000 | \n
| Alaska | \n1 | \n$17,008 | \n
| Maine | \n1 | \n$16 | \n
On our NOAA agency pulse page, we are tracking the new awards made by the agency. Unlike many of the agencies we track, NOAA does not provide up-to-date grant listings, so we use USAspending data to track new awards, which can lag by over a month. Nonetheless, we can see NOAA is far behind in its distribution of Congressionally-appropriated funds. Unlike FY 2025, when many NOAA grants were distributed in the closing days of the Biden administration, NOAA was slow to start this year and has not made new grants since mid-August. It has distributed 61% fewer funds and 74% fewer grants than the 2021\u20132024 average. We will see if the agency can catch up in the remaining weeks of the fiscal year.
\n
We will continue to track the changes and disruptions at NOAA. Keep an eye on the agency page, or follow our updates or newsletter. As with all our work at Grant Witness, crowdsourced information is key to building and verifying our databases. Please continue to submit reports of grant terminations or reinstatements using our web form, and don't hesitate to contact us.
","doi":"https://doi.org/10.59350/zbrpd-kbv50","guid":"https://doi.org/10.59350/zbrpd-kbv50","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1790553600,"rid":"cyhts-3s181","summary":"At Grant Witness, we have been tracking the damage of funding terminations across federal agencies since March 2025. We have focused on our core research agencies (NIH and NSF) and expanded to other agencies in science and public health such as EPA and CDC.","tags":["Update","National Oceanic And Atmospheric Administration"],"title":"Grant Witness is tracking NOAA Terminations and Funding","updated_at":1790621961,"url":"https://grantwitness.org/noaa/updates/2026-08-24-NOAA-targeted-ncrs","version":"v1"}},{"document":{"authors":[{"contributor_roles":[],"family":"Jivkova","given":"Katrin"}],"blog":{"authors":null,"community_id":"1e54953c-587a-4743-89af-ebd2cb45660c","created":1789516800,"current_feed_url":null,"description":null,"doi":"https://doi.org/10.59350/pensoft","favicon":"https://rogue-scholar.org/api/communities/1e54953c-587a-4743-89af-ebd2cb45660c/logo","feed_format":"application/atom+xml","feed_url":"https://blog.pensoft.net/feed/atom/","filter":null,"generator":"WordPress","home_page_url":"https://blog.pensoft.net/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"pensoft","status":"active","subfield":"2303","title":"Pensoft Blog","updated":1790604504,"use_api":true},"blog_name":"Pensoft Blog","blog_slug":"pensoft","content_html":"\n\nPensoft Publishers is pleased to announce the publication of the first issue of the Journal of Scientific and Technological Innovations (STI), an international journal for cutting-edge research across the scientific and technological spectrum.
\n\n\n\nThe journal notably adopts Pensoft's ARPHA platform, an innovative end-to-end full-featured publishing solution that streamlines the entire publication workflow \u2013 from initial manuscript submission to indexing, distribution, hosting and archiving through one coordinated platform. To date, ARPHA provides \"a home\" to about 100 open-access journals across scientific fields.
\n\n\nThe Journal of Scientific and Technological Innovations is an international, multidisciplinary academic journal supervised by Ho Chi Minh City University of Industry and Trade (HUIT). The journal publishes original, peer-reviewed research contributions.
\n\n\n\nIts editorial purview encompasses a broad range of scientific and technological domains:
\n\n\n\nThe journal is led by an international team of scholars:
\n\n\n\nEditor-in-Chief: Thi Hong Anh Le, Ho Chi Minh City University of Industry and Trade, Vietnam
\n\n\n\nDeputy Editors-in-Chief: Tuyen Truong, RMIT University, Australia; Thanh-Binh Nguyen, National Kaohsiung University of Science and Technology, Taiwan
\n\n\n\nThe first issue of STI presents six contributions spanning food science, materials science, and applied food technology:
\n\n\n\nA cost-benefit analysis of small-holder floating-house fish farming in Chau Doc City, Vietnamese Mekong Delta
\n\n\n\nA study of 35 Mekong Delta fish farmers reveals significant disparities in farm profitability across different river reaches, with operations facing demographic vulnerabilities and heavy reliance on unpaid labour.
\n\n\n\nRead more: Farzana H, Tran DD, Nguyen TD, Park E (2026) A cost-benefit analysis of small-holder floating-house fish farming in Chau Doc City, Vietnamese Mekong Delta. Journal of Scientific and Technological Innovations 1: e209236. https://doi.org/10.3897/sti.e209236
\n\n\n\nUltrasound-assisted pullulanase modification increases resistant starch content and alters the structure\u2013function properties of water caltrop starch
\n\n\n\nA combined ultrasound and enzymatic treatment increased resistant starch content in water caltrop starch from 21.72% to 28.15%, whilst enhancing its water and oil absorption capacities. The approach offers potential for developing resistant starch-enriched functional food ingredients.
\n\n\n\nRead more: Phan VM, Ngo TMP, Do TTN, Luyen TTH, Duong TT, Tran NM, Tran CH (2026) Ultrasound-assisted pullulanase modification increases resistant starch content and alters the structure\u2013function properties of water caltrop starch. Journal of Scientific and Technological Innovations 1: e209539. https://doi.org/10.3897/sti.e209539
\n\n\n\nEffects of Melatonin/\u03b2-cyclodextrin Treatment on Physicochemical and Bioactive Attributes of Cherry Tomatoes during Storage
\n\n\n\nMelatonin complexed with \u03b2-cyclodextrin reduced cherry tomato weight loss by up to 56% during refrigerated storage and proved more effective than free melatonin in preserving colour, acidity, and polyphenol content. The complex formation provides a viable strategy for extending shelf-life in harvested produce.
\n\n\n\nRead more: Tran THC, Dong TT, Nguyen PH,Le NH,Kuo C-H (2026) Effects of Melatonin/\u03b2-cyclodextrin Treatment on Physicochemical and Bioactive Attributes of Cherry Tomatoes during Storage. Journal of Scientific and Technological Innovations 1: e212043. https://doi.org/10.3897/sti.e212043
\n\n\n\nDesign, Functionalities, and Emerging Applications of Ti3C2Tx MXene-Based Composites: A Review
\n\n\n\nThis review synthesises current knowledge on titanium carbide-based MXenes, two-dimensional materials with exceptional conductivity and mechanical properties. MXene-based composites demonstrate applications in gas sensing, photocatalytic carbon dioxide reduction, and energy storage devices.
\n\n\n\nRead more: Ta QTH, Pham NN, Tran QN, Lu TMT, Noh JS (2026) Design, Functionalities, and Emerging Applications of Ti3C2Tx MXene-Based Composites: A Review. Journal of Scientific and Technological Innovations 1: e205045. https://doi.org/10.3897/sti.\u0435205045
\n\n\n\nChemical Composition, Bioactivities, and Cosmetic Applications of Coconut Oil and Its Derivatives
\n\n\n\nRecent evidence demonstrates that coconut oil, when appropriately applied, offers dermatological benefits, with antimicrobial, anti-inflammatory, and antioxidant properties. The research clarifies interactions between coconut oil and skin microflora and explains previous inconsistencies in published findings.
\n\n\n\nRead more: Dat NM, Phong DT, Long HB, Dat NT, Thy LTM (2026) Chemical Composition, Bioactivities, and Cosmetic Applications of Coconut Oil and Its Derivatives. Journal of Scientific and Technological Innovations 1: e205088. https://doi.org/10.3897/sti.e205088
\n\n\n\nRetention of bilberry anthocyanins in wood hemicelluloses coated microcapsules under the impact of humidity and packaging conditions
\n\n\n\nStorage at 8% relative humidity limited bilberry anthocyanin degradation to 17-22% over 24 weeks, with wood hemicelluloses providing comparable protection to conventional gum Arabic. The findings establish moisture control as important for developing shelf-stable natural functional powders with sustainable coating materials.
\n\n\n\nThe diversity of contributions in this first issue reflects STI’s multidisciplinary scope. Researchers interested in submitting manuscripts are invited to consult the journal's guidelines and submission portal.
\n\n\n\n***
\n\n\n\nAccess the full first issue here.
\n\n","doi":"https://doi.org/10.59350/12ke2-9s798","guid":"https://blog.pensoft.net/?p=21014","image":"https://blog.pensoft.net/wp-content/uploads/2026/09/SciComm-PR-2026-Instagram-2026-09-28T141749.018-1.png","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1790553600,"rid":"t1tqd-psc54","summary":"Tweet Pensoft Publishers is pleased to announce the publication of the first issue of the Journal of Scientific and Technological Innovations (STI), an international journal for cutting-edge research across the scientific and technological spectrum.","tags":["Academic Publishing","ARPHA Platform","Interdisciplinary Research","New Journal","Scientific Innovation"],"title":"New Journal of Scientific and Technological Innovations launches on ARPHA","updated_at":1790604910,"url":"https://blog.pensoft.net/2026/09/28/new-journal-of-scientific-and-technological-innovations-launches-on-arpha/","version":"v1"}},{"document":{"authors":[{"contributor_roles":[],"family":"Torres Munguia","given":"Juan Armando","url":"https://orcid.org/0000-0003-3432-6941"}],"blog":{"authors":null,"community_id":"ae893dd4-41e3-426b-8b48-6c9c75639805","created":1761523200,"current_feed_url":null,"description":"Exploring the art and science of working with development and humanitarian data","doi":"https://doi.org/10.59350/juan_torresmunguia","favicon":"https://rogue-scholar.org/api/communities/ae893dd4-41e3-426b-8b48-6c9c75639805/logo","feed_format":"application/rss+xml","feed_url":"https://juan-torresmunguia.netlify.app/blog/index.xml","filter":null,"generator":"Quarto","home_page_url":"https://juan-torresmunguia.netlify.app/blog/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"juan_torresmunguia","status":"active","subfield":"1804","title":"Data Art & Science Blog","updated":1783172007,"use_api":null},"blog_name":"Data Art & Science Blog","blog_slug":"juan_torresmunguia","content_html":"In these notes, I will show you how to use {ggplot} to create a waffle chart about the frequency of disease outbreaks in the world. The final output will look like this:
\n![](\"https://juan-torresmunguia.netlify.app/blog/posts/waffle-chart-disease-outbreaks//waffle-pandemics.png\")
The source of information is the global dataset of pandemic- and epidemic-prone disease outbreaks, whose data is freely available at the GitHub repository of the disease outbreaks project.
\nThe unit of analysis in the database is an outbreak, defined as the occurrence of at least one case of a specific disease in a country -or territory- during a particular year. Therefore, a country -or territory- cannot have more than one outbreak of the same disease in the same year, although it may experience outbreaks of different diseases within the same year. A country can only have multiple outbreaks of the same disease if they occur in different years.
\nThe dataset includes information on over 3,000 pandemic- and epidemic-prone disease outbreaks associated with more than 80 different infectious diseases, occurring globally from January 1996 to the present. The diseases are classified according to the International Classification of Diseases, 10th Revision (ICD-10), and the dataset contains information on the year, country, and disease for each outbreak.
\nThe dataset of pandemic- and epidemic-prone disease outbreaks is also part of the Humanitarian Data Exchange coordinated by the United Nations Office for the Coordination of Humanitarian Affairs (OCHA).
\nTo create the waffle chart, we will use the following R packages:
\n\nThe dataset is stored in an Excel. We will download it from the GitHub repository and load it into R:
\n# Define the GitHub raw content URL\nlink_data <- \"https://github.com/jatorresmunguia/disease_outbreak_news/raw/refs/heads/main/Humanitarian%20Data%20Exchange/data_2_share/disease_outbreaks_HDX.xlsx\"\n\n# Create a temporary file to store the data\ntemp <- tempfile(fileext = \".xlsx\")\n\n# Download the data from GitHub using authentication\nget_req <- GET(\n link_data,\n authenticate(Sys.getenv(\"GITHUB_PAT\"), \"\"),\n write_disk(path = temp)\n)\n\n# Read the Excel file into R\noutbreaks <- read_excel(temp)We will group the data by year and disease, and then filter the top 10 diseases from 1996 to 2024. We will also group the data by year and disease, and then filter the top 10 diseases from 2015 to 2024. Finally, we will create a factor variable for the top 10 diseases and the \"Other\" category:
\n# Top 10 diseases 1996-2024\ntop_10_icd104n <- outbreaks |>\n filter(Year < 2025) |>\n group_by(icd104n) |>\n summarise(freq = n(), .groups = \"drop\") |>\n arrange(-freq) |>\n slice_head(n = 10)\n\n# Only top-10 categories and Others\noutbreaks_year_disease_grouped <- outbreaks |>\n filter(Year < 2025 & Year > 2014) |>\n mutate(icd104n = ifelse(icd104n %in% top_10_icd104n$icd104n, icd104n, \"Other\")) |> # Group the diseases\n group_by(Year, icd104n) |>\n summarise(freq = n(), .groups = \"drop\") |> # Count the number of outbreaks\n arrange(-freq)\n\n# Create a factor variable for the top 10 diseases and the \"Other\" category\noutbreaks_year_disease_grouped <- outbreaks_year_disease_grouped |>\n mutate(icd104n = factor(icd104n,\n levels = c(\n \"Acute poliomyelitis, unspecified\", \"Classical cholera\", \"COVID-19, virus identified\",\n \"Dengue, unspecified\", \"Influenza due to identified zoonotic or pandemic influenza virus\", \"Measles\",\n \"Meningococcal meningitis\", \"Middle East respiratory syndrome coronavirus [MERS-CoV]\", \"Monkeypox\",\n \"Yellow fever, unspecified\", \"Other\"\n ),\n labels = c(\n \"Acute poliomyelitis\", \"Cholera\", \"COVID-19\", \"Dengue\", \"Influenza\",\n \"Measles\", \"Meningococcal meningitis\", \"MERS-CoV\", \"Monkeypox\",\n \"Yellow fever\", \"Other\"\n )\n ))# Add custom font\nfont_add_google(\"Roboto Condensed\", \"Roboto Condensed\") \n\nshowtext_auto()\n\n# Custom theme for the waffle chart\ntheme_waffle_chart <- function() {\n theme_minimal(\n base_family = \"Roboto Condensed\" # Base theme with custom font\n ) +\n # Custom theme settings\n theme(\n # Axis settings\n axis.title = element_blank(), # Remove axis title\n axis.line = element_blank(), # Remove axis line\n \n # Title settings\n plot.title.position = \"plot\", # Position of the title\n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 24,\n margin = margin(5, 0, 5, 0), # Top, right, bottom, left\n width = unit(1, \"npc\") # Width of the title, npc == 1 corresponds to the full width of the plot\n ),\n \n # Subtitle settings\n plot.subtitle = element_textbox(\n color = \"grey50\",\n face = \"bold\",\n size = 14,\n margin = margin(0, 0, 10, 0),\n width = unit(1, \"npc\")\n ),\n \n # Legend settings\n legend.position = \"top\",\n legend.title = element_blank(),\n legend.spacing.x = unit(0.2, \"cm\"),\n legend.key.spacing = unit(0.5, \"cm\"), # Spacing between legend keys\n legend.text = element_text(\n margin = margin(5, 2, 5, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 9\n ),\n legend.direction = \"horizontal\",\n legend.byrow = FALSE,\n \n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n size = 8,\n hjust = 0\n ),\n plot.caption.position = \"plot\",\n\n plot.background = element_rect(\n color = \"white\",\n fill = \"white\"\n ),\n plot.margin = margin(20, 40, 20, 40)\n )\n}\n\n# Title, subtitle, and caption for the waffle chart\ntitle_chart <- \"Pandemic- and epidemic-prone disease outbreaks in the world\"\nsubtitle_chart <- \"In 2024, fewer countries reported public health events linked to COVID-19, while outbreaks of viral diseases transmitted by vectors increased\"\ncaption_chart <- paste0(\"**Note:** Each square represents a country.\",\n \"<br>\", \n \"**Data:** A global dataset of pandemic- and epidemic-prone disease outbreaks (DOI: 10.1038/s41597-022-01797-2)\",\n \"<br>\", \n \"**Graphic:** Juan Torres Mungu\u00eda\"\n )\n\n# Define the colors to be used in the waffle chart\ncolors <- c(\n \"#B66DFF\", \"#7e0021\", \"#ff420e\", \"#579d1c\", \"#83caff\",\n \"#0084d1\", \"#aecf00\", \"#4b1f6f\", \"#FF86FF\", \"#ffd320\", \"grey80\"\n )ggplot(\n outbreaks_year_disease_grouped,\n aes(fill = icd104n, values = freq)\n) +\n geom_waffle(\n size = 0.25, # Size of the waffle squares\n n_rows = 10, # Number of rows in the waffle chart\n flip = TRUE, # Flip the waffle chart\n color = \"white\", # Color of the border of the waffle squares\n make_proportional = FALSE\n ) + # Do not make the waffle chart proportional\n facet_wrap(~Year, # Facet the waffle chart by year\n nrow = 1,\n strip.position = \"bottom\"\n ) + # Position of the facet strip\n scale_fill_manual(values = colors) +\n scale_x_discrete() +\n scale_y_continuous(\n labels = function(x) x * 10,\n expand = c(0, 0)\n ) + # Expand the y-axis\n coord_equal() +\n labs(\n title = title_chart,\n subtitle = subtitle_chart,\n caption = caption_chart,\n x = \"\",\n y = \"Total number of outbreaks\",\n fill = \"\"\n ) +\n guides(fill = guide_legend(nrow = 2 # Number of rows in the legend\n )) + \n theme_waffle_chart()showtext_opts(dpi = 320) # Set the resolution of the image 320 dpi is for high-quality images (\"retina\")\nggsave(\n \"waffle-pandemics.png\",\n dpi = 320,\n width = 14,\n height = 7,\n units = \"in\"\n)\nshowtext_auto(FALSE)@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {A Waffle Chart of Disease Outbreaks in the World},\n date = {2025-01-12},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/waffle-chart-disease-outbreaks/},\n langid = {en}\n}\nIn this article, I will guide you through the process of analyzing Mexico's 2024 National Survey of Victimization and Perception of Public Safety (Encuesta Nacional de Victimizaci\u00f3n y Percepci\u00f3n sobre Seguridad P\u00fablica, ENVIPE) using the {tidyverse} framework for survey analysis with the {srvyr} package. Specifically, we will calculate the crime prevalence rate per 100,000 inhabitants by state, according to the victim's sex ( Tasa de prevalencia delictiva por entidad federativa por cada cien mil habitantes, seg\u00fan sexo de la v\u00edctima).
\nOfficial tables with this indicator can be found here.
\nBy the end of this article, you will be able to use the {tidyverse} to:
\nENVIPE is an official survey that provides nationally and state-representative information on household crime victimization and perceptions of public safety in Mexico. Conducted annually since 2011 by the National Institute of Statistics and Geography (INEGI), ENVIPE collects data on the following crimes:
\nAll survey-related information, including questionnaires, methodology, glossaries, summaries, and microdata, is fully available at INEGI's website.
\nFirst, install and load the required packages:
\n\nMicrodata is available here in various formats (CSV, DBF, DTA, SAV, and RData). In this post, we will use the RData files.
\nload(archive_read(\n archive = \"https://www.inegi.org.mx/contenidos/programas/envipe/2024/microdatos/bd_envipe_2024_RData.zip\",\n file = \"BD_ENVIPE_2024.RData\"\n))Using the {archive} package, we can directly extract the required file from bd_envipe_2024_RData.zip.
This loads the following objects into the global environment:
\nTHogar: This object contains the interview result, the row of the selected respondent, the household factor, and other household characteristics.
TMod_Vic: This object includes information on the crimes suffered by the selected respondent and their household during the reference year, specifically during 2023. Information comes from the victimization module of the survey.
TPer_Vic1: This object includes information related to the perception of public security within the respondent's geographic area, antisocial behaviors in their immediate surroundings, and changes in habits due to fear of becoming a crime victim in 2023. Additionally, it presents information regarding the performance of authorities and trust in them. In other words, this table covers Sections IV and V of the main questionnaire.
TPer_Vic2: This object consolidates information related to household and personal victimization, particularly focusing on high-impact crimes experienced by the household (such as kidnapping, enforced disappearance, or homicide) or by the selected respondent. This table includes Sections VI and VII of the main questionnaire.
TSDem: This object contains the socio-demographic characteristics of household members, corresponding to Section III of the main questionnaire.
TSVivienda: This object contains the characteristics of the households captured through the questionnaire cover page, as well as the count of household residents and the number of households. In other words, it includes Sections I and II of the main questionnaire.
To estimate the crime prevalence rate, we use information from TPer_Vic2 (screening questions on victimization from the main questionnaire) and TMod_Vic (specific victimization module applied to the respondents who reported being a victim of at least one of the crimes included in the survey during the reference period). TPer_Vic2 will serve to characterize the total population and TMod_Vic to describe the victims and the associated crimes.
ENVIPE follows a probabilistic, three-stage, stratified, and clustered sampling design. The sampling unit is the household, while the unit of analysis is the population aged 18 and older in selected households. More details on the sampling design can be found here(in Spanish).
\nFrom the documentation on data structure(in Spanish), the key variables for specifying the survey design are:
\nID_PER: Person identifier
EST_DIS: Stratum
UPM_DIS: Primary sampling unit
FAC_ELE: Weighting factor for public security perception and victimization estimates.
We specify the survey design as follows:
\n# Identify victims from TMod_Vic, excluding\n# \"Wall graffiti, intentional vehicle damage, and other vandalism\" (BPCOD != \"03\")\nid_victims <- TMod_Vic |>\n filter(BPCOD != \"03\") |>\n pull(ID_PER)\n\n# Create victimization variable in TPer_Vic2\nTPer_Vic2 <- TPer_Vic2 |>\n mutate(victimization = ifelse(ID_PER %in% id_victims, 1, 0))\n\n# Convert FAC_ELE (weights) to numeric\nTPer_Vic2 <- TPer_Vic2 |>\n mutate(FAC_ELE = as.numeric(FAC_ELE))\n\n# Label SEXO variable\nTPer_Vic2 <- TPer_Vic2 |>\n mutate(SEXO = recode(SEXO, \"1\" = \"Man\", \"2\" = \"Woman\"))\n\n# Create survey object\nenvipe_design <- TPer_Vic2 |>\n as_survey_design(\n weights = FAC_ELE, # Weights (inverse probability)\n strata = EST_DIS, # Stratum\n ids = UPM_DIS, # Primary sampling unit\n nest = TRUE # Nesting within strata\n )To compute proportions from survey data, we use survey_mean() from {srvyr} package. From ENVIPE's documentation, we know that NOM_ENT represents the state, while SEXO corresponds to the respondent's sex.
envipe_design |>\n # `interact` function is used to calculate the proportion\n # over the interaction of state and sex\n group_by(interact(NOM_ENT, SEXO)) |>\n summarize(\n Proportion = survey_mean(victimization,\n level = 0.90, # Confidence level\n vartype = c(\"se\", \"ci\", \"cv\") # Variability indicators\n )\n ) |>\n ungroup() |>\n mutate(\n `Prevalence rate` = Proportion * 100000, # Rate per 100,000 people\n `Standard error` = Proportion_se * 100000,\n `Low 90% CI` = Proportion_low * 100000,\n `Upper 90% CI` = Proportion_upp * 100000,\n `Coefficient of variation %` = Proportion_cv * 100,\n State = str_to_sentence(NOM_ENT) # Capitalize only the first letter\n ) |>\n rename(Sex = SEXO) |>\n select(\n State, Sex, `Prevalence rate`, `Standard error`,\n `Low 90% CI`, `Upper 90% CI`, `Coefficient of variation %`\n ) |>\n kbl() |> # HTML table\n kable_styling() # Apply bootstrap theme to the tableThe output table looks like this:
\n| State | \nSex | \nPrevalence rate | \nStandard error | \nLow 90% CI | \nUpper 90% CI | \nCoefficient of variation % | \n
|---|---|---|---|---|---|---|
| Aguascalientes | \nMan | \n33710.81 | \n2172.7372 | \n30136.71 | \n37284.91 | \n6.445224 | \n
| Aguascalientes | \nWoman | \n32041.57 | \n1831.7081 | \n29028.45 | \n35054.68 | \n5.716662 | \n
| Baja California | \nMan | \n23043.35 | \n1604.9769 | \n20403.20 | \n25683.50 | \n6.965032 | \n
| Baja California | \nWoman | \n24480.64 | \n1536.1771 | \n21953.66 | \n27007.61 | \n6.275070 | \n
| Baja California Sur | \nMan | \n19329.94 | \n1244.1353 | \n17283.36 | \n21376.51 | \n6.436314 | \n
| Baja California Sur | \nWoman | \n23817.58 | \n1466.2963 | \n21405.56 | \n26229.60 | \n6.156361 | \n
| Campeche | \nMan | \n21832.87 | \n1455.8508 | \n19438.03 | \n24227.71 | \n6.668161 | \n
| Campeche | \nWoman | \n23213.74 | \n1381.0933 | \n20941.88 | \n25485.61 | \n5.949465 | \n
| Chiapas | \nMan | \n16602.99 | \n986.0052 | \n14981.03 | \n18224.94 | \n5.938722 | \n
| Chiapas | \nWoman | \n12074.36 | \n978.7862 | \n10464.28 | \n13684.44 | \n8.106318 | \n
| Chihuahua | \nMan | \n22303.22 | \n1234.5065 | \n20272.49 | \n24333.96 | \n5.535103 | \n
| Chihuahua | \nWoman | \n22804.31 | \n1181.0327 | \n20861.54 | \n24747.08 | \n5.178990 | \n
| Ciudad De Mexico | \nMan | \n34562.66 | \n1089.2714 | \n32770.84 | \n36354.49 | \n3.151584 | \n
| Ciudad De Mexico | \nWoman | \n30746.76 | \n990.3528 | \n29117.66 | \n32375.87 | \n3.220999 | \n
| Coahuila De Zaragoza | \nMan | \n18179.70 | \n1249.1422 | \n16124.89 | \n20234.51 | \n6.871084 | \n
| Coahuila De Zaragoza | \nWoman | \n20447.19 | \n1238.7111 | \n18409.54 | \n22484.84 | \n6.058099 | \n
| Colima | \nMan | \n21106.85 | \n1593.5970 | \n18485.43 | \n23728.28 | \n7.550140 | \n
| Colima | \nWoman | \n21295.28 | \n1351.1455 | \n19072.68 | \n23517.88 | \n6.344813 | \n
| Durango | \nMan | \n15240.74 | \n1133.1868 | \n13376.67 | \n17104.80 | \n7.435250 | \n
| Durango | \nWoman | \n17401.18 | \n1166.0910 | \n15482.99 | \n19319.37 | \n6.701219 | \n
| Guanajuato | \nMan | \n20057.38 | \n1428.0981 | \n17708.20 | \n22406.57 | \n7.120063 | \n
| Guanajuato | \nWoman | \n18930.15 | \n1287.0811 | \n16812.93 | \n21047.37 | \n6.799107 | \n
| Guerrero | \nMan | \n16981.89 | \n1494.6554 | \n14523.22 | \n19440.56 | \n8.801466 | \n
| Guerrero | \nWoman | \n14572.44 | \n1238.6234 | \n12534.93 | \n16609.94 | \n8.499769 | \n
| Hidalgo | \nMan | \n22397.82 | \n1338.0474 | \n20196.77 | \n24598.88 | \n5.974006 | \n
| Hidalgo | \nWoman | \n20215.02 | \n1157.0294 | \n18311.73 | \n22118.30 | \n5.723613 | \n
| Jalisco | \nMan | \n24390.13 | \n1548.7916 | \n21842.40 | \n26937.85 | \n6.350076 | \n
| Jalisco | \nWoman | \n25141.90 | \n1438.5170 | \n22775.58 | \n27508.23 | \n5.721592 | \n
| Mexico | \nMan | \n32906.59 | \n1632.3183 | \n30221.47 | \n35591.71 | \n4.960460 | \n
| Mexico | \nWoman | \n33023.78 | \n1483.3471 | \n30583.71 | \n35463.85 | \n4.491755 | \n
| Michoacan De Ocampo | \nMan | \n14405.17 | \n1031.1125 | \n12709.01 | \n16101.32 | \n7.157934 | \n
| Michoacan De Ocampo | \nWoman | \n15485.20 | \n1036.3320 | \n13780.46 | \n17189.94 | \n6.692401 | \n
| Morelos | \nMan | \n26995.43 | \n1517.6054 | \n24499.00 | \n29491.85 | \n5.621713 | \n
| Morelos | \nWoman | \n24869.74 | \n1688.5969 | \n22092.04 | \n27647.44 | \n6.789765 | \n
| Nayarit | \nMan | \n16075.56 | \n1516.6195 | \n13580.76 | \n18570.36 | \n9.434318 | \n
| Nayarit | \nWoman | \n18639.61 | \n1597.1842 | \n16012.28 | \n21266.94 | \n8.568766 | \n
| Nuevo Leon | \nMan | \n23602.12 | \n1261.9498 | \n21526.25 | \n25678.00 | \n5.346764 | \n
| Nuevo Leon | \nWoman | \n21440.72 | \n1211.6372 | \n19447.60 | \n23433.83 | \n5.651104 | \n
| Oaxaca | \nMan | \n12940.56 | \n1158.2366 | \n11035.29 | \n14845.83 | \n8.950434 | \n
| Oaxaca | \nWoman | \n13548.61 | \n1017.7020 | \n11874.51 | \n15222.70 | \n7.511489 | \n
| Puebla | \nMan | \n26519.64 | \n1342.0473 | \n24312.00 | \n28727.27 | \n5.060579 | \n
| Puebla | \nWoman | \n24176.45 | \n1124.4656 | \n22326.73 | \n26026.17 | \n4.651078 | \n
| Queretaro | \nMan | \n30297.50 | \n1620.6534 | \n27631.57 | \n32963.44 | \n5.349132 | \n
| Queretaro | \nWoman | \n25746.31 | \n1260.3524 | \n23673.06 | \n27819.56 | \n4.895273 | \n
| Quintana Roo | \nMan | \n22592.04 | \n1384.4250 | \n20314.70 | \n24869.39 | \n6.127932 | \n
| Quintana Roo | \nWoman | \n21775.64 | \n1237.8095 | \n19739.48 | \n23811.81 | \n5.684376 | \n
| San Luis Potosi | \nMan | \n24019.02 | \n1874.7438 | \n20935.11 | \n27102.93 | \n7.805247 | \n
| San Luis Potosi | \nWoman | \n23227.97 | \n1491.3111 | \n20774.80 | \n25681.14 | \n6.420324 | \n
| Sinaloa | \nMan | \n21181.44 | \n1122.0160 | \n19335.76 | \n23027.13 | \n5.297165 | \n
| Sinaloa | \nWoman | \n20562.22 | \n1086.2790 | \n18775.32 | \n22349.12 | \n5.282887 | \n
| Sonora | \nMan | \n26864.66 | \n1816.2677 | \n23876.94 | \n29852.37 | \n6.760808 | \n
| Sonora | \nWoman | \n25742.16 | \n1902.0264 | \n22613.37 | \n28870.95 | \n7.388759 | \n
| Tabasco | \nMan | \n25881.72 | \n1354.6212 | \n23653.40 | \n28110.04 | \n5.233891 | \n
| Tabasco | \nWoman | \n25522.45 | \n1091.8488 | \n23726.38 | \n27318.51 | \n4.277994 | \n
| Tamaulipas | \nMan | \n17788.15 | \n1075.8038 | \n16018.47 | \n19557.82 | \n6.047869 | \n
| Tamaulipas | \nWoman | \n18368.67 | \n1017.1639 | \n16695.46 | \n20041.88 | \n5.537493 | \n
| Tlaxcala | \nMan | \n23968.33 | \n1621.2198 | \n21301.47 | \n26635.20 | \n6.764007 | \n
| Tlaxcala | \nWoman | \n22295.82 | \n1473.2391 | \n19872.38 | \n24719.26 | \n6.607692 | \n
| Veracruz De Ignacio De La Llave | \nMan | \n16365.05 | \n1055.8247 | \n14628.24 | \n18101.86 | \n6.451705 | \n
| Veracruz De Ignacio De La Llave | \nWoman | \n14712.45 | \n947.9121 | \n13153.15 | \n16271.74 | \n6.442927 | \n
| Yucatan | \nMan | \n20677.82 | \n1239.0909 | \n18639.55 | \n22716.10 | \n5.992366 | \n
| Yucatan | \nWoman | \n17805.62 | \n1141.8291 | \n15927.34 | \n19683.90 | \n6.412745 | \n
| Zacatecas | \nMan | \n17328.49 | \n1346.8466 | \n15112.96 | \n19544.02 | \n7.772440 | \n
| Zacatecas | \nWoman | \n15690.25 | \n1335.9463 | \n13492.65 | \n17887.85 | \n8.514502 | \n
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Using a Tidy Approach in {R} to Analyze {Mexico's} 2024\n {ENVIPE}},\n date = {2025-01-27},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/estimating-envipe-prevalence/},\n langid = {en}\n}\nThe Standardized Precipitation Index (SPI) quantifies how precipitation deviates from the historical average over multiple time scales (e.g., from 1 to 12 months), allowing for consistent monitoring of droughts and floods across regions and time, making it a valuable tool for early warning and risk assessment systems.
\nIn this tutorial, I will walk through the process of calculating the SPI for Somalia at the subnational level using the {SCI} package in R.
\nThe rainfall dataset used in this tutorial comes from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA) and provides monthly precipitation data for Somalia's subnational districts. It is derived from CHIRPS v2 (Climate Hazards Group InfraRed Precipitation with in situ Station data), using satellite imagery combined with ground-based observations.
\nFor each 10-day period (dekad), the dataset includes the following indicators:
\nrfh: Total rainfall (mm) over 10 daysr1h: 1-month rolling rainfall sum (mm)r3h: 3-month rolling rainfall sum (mm)rfh_avg: Long-term average rainfall for the dekad (mm)r1h_avg: Long-term average for 1-month rolling rainfall (mm)r3h_avg: Long-term average for 3-month rolling rainfall (mm)rfq: Rainfall anomaly (%)r1q: Rainfall anomaly (%)r3q: 3-month rainfall anomaly (%)The data is aggregated at the subnational level and follows administrative boundaries defined by the World Food Programme (WFP), with each unit assigned a unique code (Pcode).
Before we begin, we need to install and load the necessary R packages. We will be using the {tidyverse} package for data wrangling, as well as the {SCI} package for calculating the SPI.
\nlibrary(tidyverse) \nlibrary(SCI) # For Standardized Precipitation Index calculations (climate analysis)\nlibrary(conflicted) # For resolving conflicts between function names in different packages\nlibrary(kableExtra) \n\nlibrary(ggplot2) \nlibrary(ggtext) \nlibrary(showtext) \n\nconflicts_prefer(dplyr::select) # Prefer 'select' function from dplyr package\nconflicts_prefer(dplyr::lag) # Prefer 'lag' function from dplyr package\nconflicts_prefer(dplyr::filter) # Prefer 'filter' function from dplyr packageImportant to use {conflicted} package to specify preferred functions to avoid naming conflicts between {tidyverse} and {SCI} packages.
\nData can be found here in file som-rainfall-adm2-full.csv. For Somalia, data is available from January 1981.
# Read the rainfall data for Somalia`\n# You can also download the data from the HDX platform and read it locally\nrainfall_data <- read.csv(\"https://data.humdata.org/dataset/ed6e1b4b-8094-47e6-bdf7-f6d56fa7abb9/resource/8b333d58-d69e-418c-b5e3-dd86f12eee05/download/som-rainfall-subnat-full.csv\")We preprocess the data to extract relevant information to calculate the SPI. We extract the year and month from the date column, calculate the rainfall anomaly, and group the data by year, month, and district_code.
\nrainfall_data <- rainfall_data |>\n slice(-1) |> # Remove the first row (metadata)\n filter(adm_level == 2) |>\n mutate(\n year = year(date), # Extract year from the 'date' column\n month = month(date), # Extract month from the 'date' column\n district_code = PCODE # Rename the 'ADM2_PCODE' column to 'district_code'\n ) |>\n filter(year < 2025) |> # Filter data up to 2024\n rename(\n rainfall_monthly_total = rfh, # Rename the column 'rfh' to 'rainfall_monthly_total'\n rainfall_monthly_lt_average = rfh_avg # Rename the column 'rfh_avg' to 'rainfall_monthly_lt_average'\n ) |>\n mutate(\n rainfall_monthly_total = as.numeric(rainfall_monthly_total), # Convert rainfall total to numeric\n rainfall_monthly_lt_average = as.numeric(rainfall_monthly_lt_average) # Convert long-term average to numeric\n ) |>\n select(year, month, district_code, rainfall_monthly_total, rainfall_monthly_lt_average) # Select relevant columns\n\n# Group by year, month, and district_code, and calculate mean rainfall values by district\nrainfall_data <- rainfall_data |>\n group_by(year, month, district_code) |>\n summarise(\n rainfall_monthly_total = mean(rainfall_monthly_total, na.rm = TRUE), # Calculate the mean of rainfall total\n rainfall_monthly_lt_average = mean(rainfall_monthly_lt_average, na.rm = TRUE) # Calculate the mean of long-term average\n ) |>\n ungroup() |>\n # Optionally, you can calculate rainfall anomaly as a percentage of the long-term average\n # rainfall_anomaly = 1, indicates no anomaly\n # < 1 indicates below average\n # > 1 indicates above average\n mutate(\n rainfall_anomaly = rainfall_monthly_total / rainfall_monthly_lt_average\n )rainfall_data <- rainfall_data |>\n group_by(district_code) |> # Group by district_code\n mutate(\n spi1 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1, # Start month for SPI calculation\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1, # Start month for SPI calculation\n time.scale = 1, # Time scale for SPI calculation\n distr = \"gamma\", # Distribution for fitting\n p0 = TRUE # Use initial parameter estimates\n )\n ),\n spi3 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1,\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1,\n time.scale = 3,\n distr = \"gamma\",\n p0 = TRUE\n )\n ),\n spi6 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1,\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1,\n time.scale = 6,\n distr = \"gamma\",\n p0 = TRUE\n )\n ),\n spi12 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1,\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1,\n time.scale = 12,\n distr = \"gamma\",\n p0 = TRUE\n )\n )\n ) |>\n ungroup() |>\n arrange(district_code)The output table looks like this:
\n| year | \nmonth | \ndistrict_code | \nrainfall_monthly_total | \nrainfall_monthly_lt_average | \nrainfall_anomaly | \nspi1 | \nspi3 | \nspi6 | \nspi12 | \n
|---|---|---|---|---|---|---|---|---|---|
| 1981 | \n1 | \nSO1101 | \n0.9417476 | \n1.744876 | \n0.5397218 | \n-1.1043463 | \nNA | \nNA | \nNA | \n
| 1981 | \n2 | \nSO1101 | \n2.0194174 | \n2.027616 | \n0.9959566 | \n0.1537152 | \nNA | \nNA | \nNA | \n
| 1981 | \n3 | \nSO1101 | \n32.8349507 | \n9.696980 | \n3.3861008 | \n2.0730688 | \n2.0600381 | \nNA | \nNA | \n
| 1981 | \n4 | \nSO1101 | \n18.3139157 | \n24.722222 | \n0.7407876 | \n-0.3498783 | \n0.9740066 | \nNA | \nNA | \n
| 1981 | \n5 | \nSO1101 | \n6.7928803 | \n17.754692 | \n0.3825963 | \n-1.1682128 | \n0.3482195 | \nNA | \nNA | \n
| 1981 | \n6 | \nSO1101 | \n2.7443365 | \n4.763215 | \n0.5761522 | \n-1.1812731 | \n-1.5032149 | \n0.1598925 | \nNA | \n
| 1981 | \n7 | \nSO1101 | \n17.3592230 | \n14.047789 | \n1.2357264 | \n0.7686934 | \n-0.9048239 | \n0.3881687 | \nNA | \n
| 1981 | \n8 | \nSO1101 | \n34.8381870 | \n26.758146 | \n1.3019657 | \n0.8281462 | \n0.7648482 | \n0.7461252 | \nNA | \n
| 1981 | \n9 | \nSO1101 | \n29.1682857 | \n19.003020 | \n1.5349290 | \n1.2521098 | \n1.2306838 | \n0.0193027 | \nNA | \n
| 1981 | \n10 | \nSO1101 | \n4.1974111 | \n8.226213 | \n0.5102483 | \n-0.3816705 | \n0.7653264 | \n0.1559530 | \nNA | \n
Drought events are indicated by negative SPI values, with more negative values indicating more severe drought conditions. Conversely, positive SPI values indicate wetter-than-average conditions. For demonstrating purposes, let's visualize the SPI values for Laasqoray district in Sanaag region. One efficient way to visualize SPI values is by using a box plot with SPI categories.
\nfont_add_google(\"Roboto Condensed\", \"Roboto Condensed\")\n\nshowtext_auto()\n\n# Custom theme for the chart\ntheme_boxplot_chart <- function() {\n theme_minimal(\n base_family = \"Roboto Condensed\" \n ) +\n # Custom theme settings\n theme(\n # Axis settings\n axis.title = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 14\n ),\n axis.text = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 12\n ),\n\n # Title settings\n plot.title.position = \"plot\", \n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 24,\n margin = margin(5, 0, 5, 0), \n width = unit(1, \"npc\") \n ),\n\n # Legend settings\n legend.position = \"top\",\n legend.title = element_blank(),\n legend.key.height = unit(0.5, \"cm\"),\n legend.key.width = unit(0.5, \"cm\"),\n legend.spacing.x = unit(0.2, \"cm\"),\n legend.key.spacing = unit(0.5, \"cm\"), \n legend.text = element_text(\n margin = margin(5, 0, 5, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 10\n ),\n legend.direction = \"horizontal\",\n legend.byrow = FALSE,\n\n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n face = \"italic\",\n size = 12,\n hjust = 0\n ),\n plot.background = element_rect(\n color = \"white\",\n fill = \"white\"\n ),\n plot.margin = margin(20, 40, 20, 40)\n )\n}\n\ntitle_chart <- \"1-month Standardized Precipitation Index (SPI-1) in Laasqoray, Somalia for the period 1981-2024\"\ncaption_chart <- paste0(\"**Data:** Own estimations using Somalia rainfall Indicators at Subnational Level (HDX)\",\n \"<br>\", \n \"**Graphic:** Juan Torres Mungu\u00eda\")\n\n#\nspi_colors <- data.frame(\n ymin = c(-5, -2, -1.5, -1, 0, 1, 1.5, 2),\n ymax = c(-2, -1.5, -1, 0, 1, 1.5, 2, 5),\n label = factor(\n c(\n \"Extremely dry\", \"Severely dry\", \"Moderately dry\", \"Near normal\",\n \"Near normal\", \"Moderately wet\", \"Severely wet\", \"Extremely wet\"\n ),\n levels = c(\n \"Extremely dry\", \"Severely dry\", \"Moderately dry\", \"Near normal\",\n \"Moderately wet\", \"Severely wet\", \"Extremely wet\"\n )\n ),\n fill = c(\n \"#F76D5E\", \"#FFAD72\", \"#FFE099\", \"#99FF99\",\n \"#99FF99\", \"#99EAFF\", \"#75D3FF\", \"#3D87FF\"\n )\n)\n\nrainfall_data |>\n filter(district_code == \"SO1503\") |>\n mutate(month = factor(month, levels = 1:12, labels = month.name)) |> # Convert month to factor\n ggplot(aes(x = month, y = spi1)) +\n\n # Add box plot with SPI categories\n geom_rect(\n data = spi_colors, aes(ymin = ymin, ymax = ymax, xmin = -Inf, xmax = Inf, fill = label),\n alpha = 0.2, inherit.aes = FALSE\n ) +\n geom_boxplot(size = 1, outlier.size = 2, outlier.shape = 16) +\n scale_fill_manual(values = setNames(spi_colors$fill, spi_colors$label)) +\n labs(\n title = title_chart,\n caption = caption_chart,\n x = \"\",\n y = \"SPI-1\",\n fill = \"\"\n ) +\n guides(fill = guide_legend(nrow = 1)) + # Number of rows in the legend\n theme_boxplot_chart()showtext_opts(dpi = 320) # Set the resolution of the image 320 dpi is for high-quality images (\"retina\")\nggsave(\n \"spi-somalia.png\",\n dpi = 320,\n width = 12,\n height = 9,\n units = \"in\"\n)\nshowtext_auto(FALSE) # Turn off the showtext functionality![](\"https://juan-torresmunguia.netlify.app/blog/posts/spi_somalia//spi-somalia.png\")
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {How to Compute and Visualize the {Standardized}\n {Precipitation} {Index} {(SPI)} for Drought Analysis},\n date = {2025-02-08},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/spi_somalia/},\n langid = {en}\n}\nPreviously, in Part 1 of this tutorial, we calculated the Standardized Precipitation Index (SPI) for different time scales using the {SCI} package in R. Now, we will map the SPI values for the districts of Somalia using {sf} and {ggplot2} packages in R.
\nThe information of the Somalian sub national administrative boundaries is obtained from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA). The information comes as a shapefile .shp and is available here in a zip file named som_adm_ocha_20250108_AB_SHP.zip. The shapefile contains the administrative boundaries of Somalia at the subnational level. In particular, we will use file som_admbnda_ADM2_ocha_20250108.shp, that contains the shapefile of Somalia at the district level. I downloaded this file and located it in /shp folder.
Before we begin, we need to install and load the necessary R packages.
\n\nData can be found here in file som-rainfall-adm2-full.csv. For Somalia, data is available from January 1981.
I will add the SPI values previously calculated to the shapefile data frame. I saved the SPI values in a RData file named rainfall_data.RData.
load(\"rainfall_data.RData\")\n\n# First, we need to filter the SPI values for the last available month, i.e. December 2024\n# to have only one observation per district\nrainfall_data_dec2024 <- rainfall_data |>\n filter(month == 2 & year == 2024)\n\n# Merge the shapefile with the SPI values\nshp_somalia <- shp_somalia |>\n left_join(rainfall_data_dec2024, by = c(\"district_code\"))font_add_google(\"Roboto Condensed\", \"Roboto Condensed\")\n\nshowtext_auto()\n\ntheme_spi_map <- function() {\n theme_minimal(\n base_family = \"Roboto Condensed\" \n ) +\n # Custom theme settings\n theme(\n\n # Title settings\n plot.title.position = \"plot\", \n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 22,\n margin = margin(5, 0, 5, 0), \n width = unit(1, \"npc\") \n ),\n plot.subtitle = element_textbox(\n color = \"grey50\",\n face = \"bold\",\n size = 14,\n margin = margin(20, 0, 10, 0),\n width = unit(1, \"npc\")\n ),\n\n legend.position = \"top\",\n legend.title = element_blank(),\n legend.key.height = unit(0.1, \"cm\"), \n legend.key.width = unit(0.1, \"cm\"), \n legend.spacing.x = unit(0.1, \"cm\"),\n legend.key.spacing = unit(0.1, \"cm\"), \n legend.text = element_text(\n margin = margin(5, 0, 5, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 10\n ),\n legend.direction = \"horizontal\",\n legend.byrow = FALSE,\n\n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n face = \"italic\",\n size = 10,\n hjust = 0\n ),\n plot.background = element_rect(\n color = \"white\",\n fill = \"white\"\n ),\n plot.margin = margin(20, 40, 20, 40)\n )\n}\n\n# Title, subtitle, and caption for the waffle chart\ntitle_chart <- \"Standardized Precipitation Index (SPI) in Somalia\"\nsubtitle_chart <- \"February 2024\"\ncaption_chart <- \"**Data:** Somalia: Own estimations using Somalia Rainfall Indicators at Subnational Level (HDX) <br> **Graphic:** Juan Torres Mungu\u00eda\"Now, lets create a histogram of the SPI values for the districts of Somalia to use as a legend for the map.
\nspi_colors <- data.frame(\n ymin = c(-5, -2, -1.5, -1, 0, 1, 1.5, 2),\n ymax = c(-2, -1.5, -1, 0, 1, 1.5, 2, 5),\n label = factor(\n c(\n \"Extreme dry\", \"Severe dry\", \"Moderate dry\", \"Near normal\",\n \"Near normal\", \"Moderate wet\", \"Severe wet\", \"Extreme wet\"\n ),\n levels = c(\n \"Extreme dry\", \"Severe dry\", \"Moderate dry\", \"Near normal\",\n \"Moderate wet\", \"Severe wet\", \"Extreme wet\"\n )\n ),\n fill = c(\n \"#F76D5E\", \"#FFAD72\", \"#FFE099\", \"#E5FFE5\",\n \"#E5FFE5\", \"#99EAFF\", \"#75D3FF\", \"#3D87FF\"\n )\n)\n\nshp_somalia <- shp_somalia %>%\n mutate(spi_category = cut(spi1,\n breaks = c(-Inf, spi_colors$ymax),\n labels = spi_colors$label,\n right = TRUE\n )) |> # right = TRUE means that the intervals are closed on the right\n filter(!is.na(spi1))\n\nhistogram_spi <- shp_somalia |>\n filter(!is.na(spi_category)) |>\n ggplot(aes(x = spi_category, fill = spi_category)) +\n geom_bar(color = \"white\") +\n scale_fill_manual(values = setNames(spi_colors$fill, spi_colors$label)) +\n labs(\n title = \"\",\n x = \"\",\n y = \"Number of districts\"\n ) +\n theme_spi_map() +\n theme(\n legend.position = \"none\",\n axis.title = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 7\n ),\n axis.text = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 7\n ),\n axis.line = element_blank(),\n panel.grid = element_blank()\n ) +\n coord_flip()# Read the shapefile of Somalia with the interior region borders\nshp_somalia_region <- read_sf(\"shp/som_admbnda_adm1_ocha_20250108.shp\")\n\nshp_somalia_region <- st_transform(shp_somalia_region, 3857) # Web Mercator projection\nshp_somalia <- st_transform(shp_somalia, 3857)\n\nmap_spi <- ggplot() +\n geom_sf(\n data = shp_somalia, aes(fill = spi_category), color = \"grey\", linewidth = 0.5\n ) +\n scale_fill_manual(values = setNames(spi_colors$fill, spi_colors$label)) +\n geom_sf(data = shp_somalia_region, fill = NA, color = \"black\", linewidth = 0.9) +\n labs(\n title = title_chart,\n subtitle = subtitle_chart,\n caption = caption_chart,\n x = \"\",\n y = \"\",\n fill = \"\"\n ) +\n theme_spi_map() +\n theme(\n legend.position = \"none\",\n # Axis settings\n axis.title = element_blank(),\n axis.line = element_blank(),\n axis.text = element_blank(),\n axis.ticks = element_blank(),\n panel.grid = element_blank()\n )map_spi +\n inset_element(histogram_spi, left = 0.4, bottom = 0, right = 1, top = 0.4, on_top = FALSE)showtext_opts(dpi = 320) \nggsave(\n \"spi-somalia-map.png\",\n dpi = 320,\n width = 12,\n height = 9,\n units = \"in\"\n)\nshowtext_auto(FALSE)![](\"https://juan-torresmunguia.netlify.app/blog/posts/spi_somalia_visualization//spi-somalia-map.png\")
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {How to Compute and Visualize the {Standardized}\n {Precipitation} {Index} {(SPI)} for Drought Analysis},\n date = {2025-02-14},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/spi_somalia_visualization/},\n langid = {en}\n}\nACLED (Armed Conflict Location & Event Data Project) provides detailed data on political violence and protest events across the globe. This blog post explores the visualization of armed conflict in the Lake Chad Basin, focusing on the dynamics of violence against civilians in the region. Using ACLED's data in combination with shapefiles from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA), we leverage {ggplot2}, {sf} and the {tidyverse} to create insightful maps that highlight the impact of armed conflict.
\nLoad necessary libraries for data processing and visualization.
\nlibrary(tidyverse) # For data manipulation and visualization\nlibrary(sf) # For handling spatial data\nlibrary(ggplot2) # For data visualization\nlibrary(patchwork) # For combining plots\nlibrary(ggtext) # Use rich text elements in ggplot\nlibrary(showtext) # Add custom fonts to plots\nlibrary(ggspatial) # For interacting with spatial data using ggplot2\nlibrary(rnaturalearth) # To provide map data that can be visualized using other R packages\nlibrary(rnaturalearthdata) # To store data for the rnaturalearth package\n\n# Set environment to English (useful for date/time formatting)\ninvisible(Sys.setenv(LANG = \"en\"))\ninvisible(Sys.setlocale(\"LC_TIME\", \"English\"))Accessing ACLED data requires registering on their ACLED Access Portal to generate a unique access key. Once registered, you can download the data using their Data Export Tool. I downloaded all the events occurred in Africa since 1997 and save them in file acled_africa_1997-2025.csv. To have a better understanding of their data, they provide a codebook.
conflict_data <- read.csv(\"acled_africa_1997-2025.csv\")Information of the national and sub national administrative boundaries in the Lake Chad Basin is obtained from the HDX. The information comes as a shapefile .shp and is available here in a zip file named lcb_admbnda_adm2_ocha.zip. The shapefile contains the administrative boundaries of Lake Chad Basin region at the subnational level. In particular, we will use file lcb_admbnda_adm2_ocha.shp. I also downloaded the administrative boundaries of the West and Central Africa region available here.
# Lake Chad Basin Area - Administrative boundaries for Nigeria, Chad, Cameroon and Niger\nshp_lake_chad_basin <- read_sf(\"lcb_admbnda_adm2_ocha/lcb_admbnda_adm2_ocha.shp\")\n\n# West and Central Africa - Administrative boundaries levels 0 - 2 and Settlements\nwest_central_africa <- read_sf(\"wca_admbnda_adm2_ocha/wca_admbnda_adm2_ocha.shp\")As indicated in the ACLED's codebook, geographic location of the events is specified in columns latitude and longitude in the conflict_data object and they are using the EPSG:4326 coordinate reference system. We have to ensure, that both shp_lake_chad_basin and west_central_africa have the same coordinate reference system. To achieve this, we use st_transform() from {sf} package.
# st_transform to convert coordinates to ACLED's data coordinate reference system\nshp_lake_chad_basin <- st_transform(shp_lake_chad_basin,\n crs = 4326\n)\nwest_central_africa <- st_transform(west_central_africa,\n crs = 4326\n)Now, we subset the events to exclusively account for those occurred between 2020 and 2024 in the Lake Chad Basin. To achieve this, we use the {tidyverse} package for data wrangling and {sf} for spatial filtering. By applying the st_filter function, we retain only events located within the Lake Chad Basin area.
conflict_data <- conflict_data |>\n # Keep only events from 2020 to 2024\n filter(year >= 2020 & year < 2025) |>\n # Select only \"Violence against civilians\" events of type \"Attack\"\n filter(event_type == \"Violence against civilians\" & sub_event_type == \"Attack\") |>\n # Keep events from Nigeria, Niger, Cameroon, and Chad\n filter(country %in% c(\"Nigeria\", \"Niger\", \"Cameroon\", \"Chad\")) |>\n # Extract the month name from event_date\n mutate(\n month = month(dmy(event_date),\n label = TRUE\n )\n ) |>\n # Select only relevant columns\n select(\n event_id_cnty, event_date, month, year, event_type, sub_event_type,\n actor1, actor2, location, latitude, longitude, fatalities, notes\n )\n\n# Convert the data to a spatial object using latitude and longitude (WGS 84 CRS)\nconflict_data <- st_as_sf(conflict_data,\n coords = c(\"longitude\", \"latitude\"),\n crs = 4326\n)\n\n# Filter events that fall within the Lake Chad Basin shapefile\nconflict_data_lake_chad <- conflict_data |>\n st_filter(shp_lake_chad_basin)We then processes administrative boundary data to extract relevant regions in West and Central Africa, focusing on the Lake Chad Basin area.
\n# Filter to include only countries in the Lake Chad Basin region\n# I included Central African Republic for visualization purposes\nshp_countries <- west_central_africa |>\n filter(admin0Name %in% c(\"Nigeria\", \"Niger\", \"Chad\", \"Cameroon\", \"Central African Republic\"))\n\n# Create country codes and format administrative division codes\nshp_countries <- shp_countries |>\n mutate(\n country_code = case_when(\n admin0Name == \"Nigeria\" ~ \"NGA\",\n admin0Name == \"Niger\" ~ \"NER\",\n admin0Name == \"Chad\" ~ \"TCD\",\n admin0Name == \"Cameroon\" ~ \"CMR\"\n ),\n # Extract and format adm1 and adm2 codes\n adm1_code = str_pad(str_extract(admin1Pcod, \"\\\\d+\"), 3, pad = \"0\"),\n adm2_code = str_pad(str_sub(admin2Pcod, 5, 6), 3, pad = \"0\"),\n # Create a unique administrative code by combining country and division codes\n adm_code = paste0(country_code, adm1_code, adm2_code)\n ) |>\n # Rename columns for clarity\n rename(\n country = admin0Name,\n adm1 = admin1Name,\n adm2 = admin2Name\n ) |>\n # Select only relevant columns\n select(adm_code, country, adm1, adm2, geometry)\n\n# Identify administrative units within the Lake Chad Basin\nshp_countries <- shp_countries |>\n mutate(lake_chad_basin = ifelse(adm_code %in% shp_lake_chad_basin$Rowcacode2,\n \"Yes\", \"No\"\n ))\n\n# Subset to retain only the administrative units inside the Lake Chad Basin\nshp_lake_chad <- shp_countries |>\n filter(lake_chad_basin == \"Yes\")\n\n# Compute the bounding box for visualization\nbbox <- st_bbox(shp_lake_chad)\nx_min <- bbox[\"xmin\"]\nx_max <- bbox[\"xmax\"]\ny_min <- bbox[\"ymin\"]\ny_max <- bbox[\"ymax\"]\n\n# Aggregate administrative level 1 (regions/states) within the Lake Chad Basin\nshp_adm1_aggregated <- shp_countries |>\n group_by(adm1, lake_chad_basin) |>\n summarise(geometry = st_union(geometry))\n\n# Keep only regions in the Lake Chad Basin and clean names\nshp_adm1_aggregated <- shp_adm1_aggregated |>\n filter(lake_chad_basin == \"Yes\") |>\n mutate(adm1 = str_replace_all(adm1, \"_\", \" \")) |>\n # Compute centroid coordinates for labeling\n mutate(\n x = st_coordinates(st_centroid(geometry))[, 1],\n y = st_coordinates(st_centroid(geometry))[, 2]\n )Now, we create a hexagonal grid over the Lake Chad region with a cell size of approximately 10.69 km and fatalities and attack counts are aggregated per grid cell.
\n# Reproject Lake Chad Basin shapefile to UTM Zone 33N (EPSG:32633)\n# This is important to set a cell size in km\nshp_lake_chad <- st_transform(shp_lake_chad, crs = 32633)\n\n# Create a hexagonal grid over the Lake Chad Basin with a cell size of ~10.69 km\nnc_grid <- st_make_grid(shp_lake_chad,\n cellsize = 10690,\n square = FALSE\n)\n\n# Filter the grid to keep only hexagons intersecting the Lake Chad Basin\nnc_grid_filter <- nc_grid |>\n st_sf() |>\n st_filter(shp_lake_chad)\n\n# Reproject the grid and Lake Chad Basin back to WGS84 (EPSG:4326)\nnc_grid_filter <- st_transform(nc_grid_filter, crs = 4326)\nshp_lake_chad <- st_transform(shp_lake_chad, crs = 4326)\n\n# Spatially join conflict data with the hexagonal grid\nconflict_data_grid <- st_join(nc_grid_filter, conflict_data_lake_chad, left = FALSE)\n\n# Aggregate conflict data by grid cell, summing fatalities and counting attacks\nfatalities_per_grid <- conflict_data_grid |>\n group_by(nc_grid) |>\n summarise(\n fatalities = sum(fatalities, na.rm = TRUE),\n num_attacks = n()\n ) |>\n st_as_sf() |>\n # Compute centroid coordinates for visualization\n mutate(\n x = st_coordinates(st_centroid(nc_grid))[, 1],\n y = st_coordinates(st_centroid(nc_grid))[, 2]\n )We need to relocate the coordinates of the labels for specific regions for better readability.
\n# Adjust label positions for specific administrative regions to improve readability\nshp_adm1_aggregated <- shp_adm1_aggregated |>\n mutate(\n x = case_when(\n adm1 == \"Extreme Nord\" ~ x - 0.1, # Shift label slightly left\n TRUE ~ x\n ),\n y = case_when(\n adm1 == \"Extreme Nord\" ~ y - 0.5, # Move label downward\n adm1 == \"Adamawa\" ~ y + 0.3, # Move label upward\n adm1 == \"Borno\" ~ y - 0.4, # Adjust label downward\n TRUE ~ y\n )\n )\n\n# Aggregate country geometries by merging administrative boundaries\n# Dissolve some shapes\nshp_countries_countries <- shp_countries |>\n group_by(country) |>\n summarise(geometry = st_union(geometry)) |>\n st_cast(\"POLYGON\") |>\n group_by(country) |>\n # Retain the largest polygon for each country\n filter(st_area(geometry) == max(st_area(geometry))) \n\n# Define manual label positions for country names\ncountry_labels <- tibble(\n country = c(\"Nigeria\", \"Niger\", \"Chad\", \"Cameroon\"),\n x = c(10.15, 11, 15, 13.6), # X-coordinates for label placement\n y = c(10, 15, 14.2, 9) # Y-coordinates for label placement\n)We first create an empty map of the countries and highlight the Lake Chad Basin area using {ggplot2}.
\n# Automatically use showtext for handling fonts\nshowtext_auto()\n\n# Add Google font \"Voltaire\"\nfont_add_google(\"Voltaire\", \"Voltaire\")\n\nlake_chad_map <- ggplot() +\n\n # Plot country borders with a light fill color\n geom_sf(\n data = shp_countries,\n fill = \"#ECEFF1\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n ) +\n\n # Plot Lake Chad Basin area with a darker fill color\n geom_sf(\n data = shp_lake_chad,\n fill = \"#E0E0E0\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n )Now, we add the fatalities information.
\nlake_chad_map <- lake_chad_map +\n # Plot fatalities per grid with a color gradient\n geom_sf(\n data = fatalities_per_grid,\n aes(fill = fatalities),\n color = NA,\n alpha = 1\n ) +\n scale_fill_gradient(name = \"Fatalities\", \n low = \"#FCC5C0\", high = \"#67001F\")We add different color to subnational and national boundaries and add their corresponding labels.
\nlake_chad_map <- lake_chad_map +\n # Plot aggregated administrative region boundaries with white lines\n geom_sf(\n data = shp_adm1_aggregated,\n fill = NA,\n color = \"white\",\n linewidth = 0.5\n ) +\n\n # Plot country boundaries with black lines\n geom_sf(\n data = shp_countries_countries,\n fill = NA,\n color = \"black\",\n linewidth = 1.25\n ) +\n # Add labels for administrative regions (e.g., \"Extreme Nord\")\n geom_text(\n data = shp_adm1_aggregated,\n aes(\n x = x, y = y,\n # Wrap text if needed for long region names\n label = str_wrap(adm1, 1) \n ),\n colour = \"#999\",\n fontface = \"bold\",\n size = 3,\n family = \"Voltaire\"\n ) +\n \n # Add labels for countries with a larger font\n geom_label(\n data = country_labels,\n aes(x = x, y = y, label = country),\n colour = \"black\",\n fontface = \"bold\",\n size = 4,\n family = \"Voltaire\"\n )Finally, we center the map at the Lake Chad Basin area and add format for a better visualization.
\n# Title, subtitle, and caption for the waffle chart\ntitle_chart <- \"Armed conflict in the Lake Chad Basin, 2020-2024\"\nsubtitle_chart <- \"Fatalities in attacks perpetrated by organized armed actors against civilians\"\n\nlake_chad_map <- lake_chad_map +\n\n labs(\n title = title_chart,\n subtitle = subtitle_chart,\n caption = \"**Data:** Armed Conflict Location & Event Data (ACLED) <br> **Graphic:** Juan Torres Mungu\u00eda\",\n x = \"\",\n y = \"\",\n fill = \"\"\n ) +\n\n # Set the map's coordinate limits and prevent expansion\n coord_sf(\n xlim = c(x_min, x_max),\n ylim = c(y_min, y_max),\n expand = FALSE\n ) +\n\n # Add a scale bar in the bottom left\n annotation_scale(location = \"bl\", width_hint = 0.2) +\n\n # Add a north arrow in the top left corner\n annotation_north_arrow(\n location = \"tl\",\n which_north = \"true\",\n style = north_arrow_fancy_orienteering,\n height = unit(1, \"cm\"),\n width = unit(1, \"cm\")\n ) +\n\n # Customize the legend\n guides(\n fill = guide_legend(position = \"bottom\")\n ) +\n\n # Set the base theme without any background elements\n theme_void(base_family = \"Voltaire\") +\n\n # Customize the map's text and legend styling\n theme(\n # Title settings\n plot.title.position = \"plot\", # Position of the title\n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 22,\n margin = margin(5, 3, 5, 3), # Top, right, bottom, left\n # Width of the title, npc == 1 corresponds to the full width of the plot\n width = unit(1, \"npc\") \n ),\n plot.subtitle = element_textbox(\n color = \"grey50\",\n face = \"bold\",\n size = 14,\n margin = margin(20, 3, 10, 3),\n width = unit(1, \"npc\")\n ),\n\n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n size = 10,\n margin = margin(10, 3, 5, 3)\n ),\n legend.title = element_text(hjust = 0.5, color = \"black\", face = \"bold\"),\n legend.text = element_text(color = \"black\", face = \"bold\"),\n legend.direction = \"horizontal\",\n legend.title.position = \"top\",\n legend.text.position = \"bottom\",\n legend.key.size = unit(0.65, units = \"cm\"),\n legend.background = element_rect(fill = NA, color = NA),\n legend.margin = margin(t = 5, r = 10, b = 5, l = 10),\n plot.background = element_rect(fill = \"#ECEFF1\", color = NA)\n )We create a zoomed-in map of the Lake Chad Basin inset within the African continent, using ggplot2. It includes the ne_countries function from the {rnaturalearth} package to retrieve the map of Africa. The map of the Lake Chad Basin is overlaid on top of Africa, and a rectangle is drawn to highlight the region of interest (Lake Chad Basin). The map is customized with a white background and black borders.
\n# Fetch African countries map from the rnaturalearth package\nafrica <- ne_countries(continent = \"Africa\", returnclass = \"sf\")\n\n# Create a zoomed-in map of Lake Chad Basin within the African map\nlake_chad_map_zoom <- ggplot() +\n\n # Plot the Africa map with a light fill color\n geom_sf(\n data = africa,\n fill = \"#ECEFF1\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n ) +\n\n # Plot Lake Chad Basin area with a darker fill color\n geom_sf(\n data = shp_lake_chad,\n fill = \"#E0E0E0\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n ) +\n\n # Add a rectangle highlighting the Lake Chad Basin on the Africa map\n geom_rect(\n aes(xmin = x_min, xmax = x_max, ymin = y_min, ymax = y_max),\n color = \"black\",\n fill = NA,\n linewidth = 0.2\n ) +\n\n # Remove grid and axes for a clean map appearance\n theme_void() +\n\n # Customize borders and background of the map\n theme(\n panel.border = element_rect(\n color = \"black\",\n fill = NA,\n size = 0.8\n ),\n plot.background = element_rect(fill = \"white\")\n )Finally, we combine the main map of the Lake Chad Basin with an inset map of Africa showing the location of the Lake Chad Basin. The inset map is positioned within the main map using the inset_element() function of the {patchwork} package. Additionally, we set the resolution of the saved image to 320 dpi (high resolution) for better quality. Finally, the map is saved as a PNG file using ggsave()
# locate main and inset map\nlake_chad_map +\n inset_element(\n align_to = \"plot\",\n lake_chad_map_zoom,\n right = 1,\n left = 0.84,\n top = 1,\n bottom = 0.91,\n on_top = TRUE\n )# Set the resolution of the image 320 dpi is for high-quality images (\"retina\")\nshowtext_opts(dpi = 320) \nggsave(\n \"lake-chad-map.png\",\n dpi = 320,\n width = 5,\n height = 12,\n units = \"in\"\n)\nshowtext_auto(FALSE)![](\"https://juan-torresmunguia.netlify.app/blog/posts/visualizing-conflict-lake-chad-basin/lake-chad-map.png\")
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Visualizing Armed Conflict Using {ACLED's} Data and {HDX}\n Shapefiles},\n date = {2025-03-09},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/visualizing-conflict-lake-chad-basin/},\n langid = {en}\n}\nWe begin by loading the necessary R packages:
\n\nWe retrieve the latest dataset of disease outbreaks from GitHub:
\n# Define the API URL from GitHub containing the latest dataset update\nurl_api <- \"https://api.github.com/repos/jatorresmunguia/disease_outbreak_news/contents/Last%20update\"\n\n# Retrieve the list of files in the \"Last update\" folder\nfile_list <- fromJSON(content(GET(url_api), as = \"text\"))$name\n\n# Filter the most recent file that starts with \"outbreaks\"\nlast_file <- file_list[grepl(file_list, pattern = paste0(\"^outbreaks\"))]\n\n# Select the file with a .csv extension\nrdata_file <- last_file[grepl(\".csv$\", last_file)]\n\n# Construct the CSV file URL from GitHub and read the data into a dataframe\nlast_version <- read.csv(paste0(\"https://raw.githubusercontent.com/jatorresmunguia/disease_outbreak_news/refs/heads/main/Last%20update/\", rdata_file), header = TRUE)\n\n# Remove unnecessary ICD-11 classification columns\nlast_version <- last_version |> select(-c(icd11c1, icd11c2, icd11c3, icd11l1, icd11l2, icd11l3))We integrate additional geographic classifications from WHO and UNSD.
\n# Define the URL for the WHO regional classification dataset\nurl_who_class <- \"https://cdn.who.int/media/docs/default-source/air-pollution-documents/air-quality-and-health/un-agencies-region-classification-for-country.xlsx?sfvrsn=289af35f_3\"\n\n# Create a temporary file to store the downloaded Excel file\ntemp_file <- tempfile(fileext = \".xlsx\")\n\n# Download the WHO classification file and save it to the temporary file\nGET(url_who_class, write_disk(temp_file, overwrite = TRUE))Response [https://cdn.who.int/media/docs/default-source/air-pollution-documents/air-quality-and-health/un-agencies-region-classification-for-country.xlsx?sfvrsn=289af35f_3]\n Date: 2025-04-28 06:21\n Status: 200\n Content-Type: application/vnd.openxmlformats-officedocument.spreadsheetml.sheet\n Size: 28.9 kB\n<ON DISK> C:\\Users\\jator\\AppData\\Local\\Temp\\RtmpiGD1DX\\file990423923f25.xlsx# Read the WHO classification data from the downloaded Excel file\nwho_class <- read_excel(temp_file)\n\n# Select relevant columns (ISO country code and WHO region name)\n# Rename columns for consistency with the main dataset\nlast_version <- who_class |>\n select(`ISO Country code`, `WHO Region name2`) |>\n rename(iso3 = `ISO Country code`, who_region = `WHO Region name2`) |>\n right_join(last_version, by = \"iso3\") # Merge with the main dataset based on the ISO3 country code# Read the UNSD regional classification dataset from an Excel file\n# The dataset is available at: https://unstats.un.org/unsd/methodology/m49/\nunsd_class <- read_excel(path = \"UNSD \u2014 Methodology.xlsx\")\n\n# Select relevant columns: ISO-alpha3 country code, region name, and sub-region name\n# Rename columns to match the main dataset\nlast_version <- unsd_class |>\n select(`ISO-alpha3 Code`, `Region Name`, `Sub-region Name`) |>\n rename(\n iso3 = `ISO-alpha3 Code`,\n unsd_region = `Region Name`,\n unsd_subregion = `Sub-region Name`\n ) |>\n right_join(last_version, by = \"iso3\") # Merge with the main dataset based on the ISO3 country codeTo facilitate better organization, we:
\n# Arrange the dataset in descending order by Year, then by ISO3 country code, and ICD-10 4-character code\nlast_version <- last_version |>\n arrange(desc(Year), iso3, icd104c) |>\n # Create a unique outbreak ID by concatenating Year, ISO3 country code, and ICD-10 4-character code\n mutate(id_outbreak = paste0(Year, iso3, icd104c)) |>\n # Add an empty row at the beginning of the dataset (for later HXL hashtag insertion)\n add_row(.before = 1)We define HXL hashtags for each relevant column:
\n# Define the mapping of column names to their corresponding HXL hashtags\ncols2hxl <- c(\n id_outbreak = \"#Year+iso3+icd4\", # Unique outbreak ID\n Country = \"#country+name\", # Country name\n iso2 = \"#country+code+iso2\", # ISO 2-letter country code\n iso3 = \"#country+code+iso3\", # ISO 3-letter country code\n unsd_region = \"#region+unsd\", # UNSD region\n unsd_subregion = \"#subregion+unsd\", # UNSD subregion\n who_region = \"#region+who\", # WHO region\n Year = \"#date+year\", # Year of the outbreak\n icd10n = \"#disease+name+icd\", # Disease name (ICD-10)\n icd103n = \"#disease+name+icd3\", # Disease name (ICD-10 3-character)\n icd104n = \"#disease+name+icd4\", # Disease name (ICD-10 4-character)\n icd10c = \"#disease+code+icd\", # Disease code (ICD-10)\n icd103c = \"#disease+code+icd3\", # Disease code (ICD-10 3-character)\n icd104c = \"#disease+code+icd4\", # Disease code (ICD-10 4-character)\n Disease = \"#disease+name\", # Disease name\n DONs = \"#news+id\", # News ID\n Definition = \"#x_disease+definition\" # Disease definition\n)\n\n# Add HXL hashtags to the relevant columns in the dataset\nlast_version <- last_version |>\n mutate(across(all_of(names(cols2hxl)), ~ replace(., 1, cols2hxl[cur_column()]))) |>\n # Select the relevant columns for the final dataset\n select(id_outbreak, Year, icd10n, icd103n, icd104n, icd10c, icd103c, icd104c, Disease, Definition, Country, iso2, iso3, unsd_region, unsd_subregion, who_region, DONs)Finally, we save the dataset as an Excel file:
\nwrite_xlsx(x = list(Data = last_version), path = \"disease_outbreaks_HDX.xlsx\")@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Working with {Humanitarian} {Exchange} {Language} {(HXL)}\n {Hashtags}},\n date = {2025-03-26},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/hdx-tags-to-outbreaks-data/},\n langid = {en}\n}\nAccording to the 2021-22 Domestic Violence module of the Kingdom of Cambodia's Demographic and Health Survey (CDHS), over 21% of Cambodian women aged 15\u201349 who have ever had a husband (or intimate partner), have experienced physical, sexual, economic, and/or emotional intimate partner violence in their lifetime. This figures varies significantly across the country, with some provinces reporting rates as high as 45% while others report rates as low as 6%.
\nA better understanding of the geographical distribution of intimate partner violence (IPV) in Cambodia can help policymakers and organizations to target their efforts more effectively. In this blog post, I will show how to create choropleth maps using R and shapefiles to visualize the prevalence of IPV in Cambodia at the provincial level. Additionally, I will use the {patchwork} package to combine these maps into the same graphic. The result will look like this:
\n![](\"https://juan-torresmunguia.netlify.app/blog/posts/map-intimate-partner-violence-cambodia//ipv-cambodia-map.png\")
Information comes from the document \"Further Data Analysis Report. From the Cambodia Demographic and Health Survey 2021-2022\", available here. In particular, I will use Table A.1: Intimate-partner violence by background characteristics.
\nFirst, we install and load the necessary R packages.
\nlibrary(tidyverse) \nlibrary(sf) # For handling shapefiles\nlibrary(readxl) # To read Excel files\nlibrary(kableExtra) # For data tables\nlibrary(ggplot2) \nlibrary(ggtext) \nlibrary(showtext) \nlibrary(patchwork) # To combine multiple ggplots\nlibrary(colorspace) # For color palettes\nlibrary(grid) # To add lines and text elements to graphical objectsI saved the information from Table A.1: Intimate-partner violence by background characteristics in the Excel file intimate-partner-violence-cambodia.xlsx.
\nipv_cambodia <- read_excel(\"intimate-partner-violence-cambodia.xlsx\")This information looks like this:
\nipv_cambodia |>\n kbl(\n caption = \"Lifetime experience of specific acts of intimate partner physical violence committed by their current or most recent partner by provinces, Cambodia 2021-22\"\n ) |>\n kable_paper(\"hover\", full_width = F)| Region | \nEmotional violence | \nPhysical violence | \nSexual violence | \nPhysical or sexual | \nPhysical or sexual or emotional | \nNumber of women who ever had a husband/intimate partner | \n
|---|---|---|---|---|---|---|
| Banteay Meanchey | \n31.3 | \n18.2 | \n4.6 | \n20.2 | \n34.9 | \n235 | \n
| Battambang | \n24.9 | \n16.1 | \n6.3 | \n19.8 | \n30.4 | \n378 | \n
| Kampong Cham | \n7.7 | \n5.2 | \n1.1 | \n5.8 | \n11.5 | \n352 | \n
| Kampong Chhnang | \n19.2 | \n6.7 | \n1.7 | \n6.7 | \n19.2 | \n213 | \n
| Kampong Speu | \n16.6 | \n9.5 | \n2.6 | \n10.1 | \n19.1 | \n362 | \n
| Kampong Thom | \n32.4 | \n12.1 | \n5.0 | \n14.4 | \n35.1 | \n242 | \n
| Kampot | \n20.8 | \n6.4 | \n4.0 | \n9.0 | \n23.0 | \n218 | \n
| Kandal | \n23.9 | \n6.2 | \n2.0 | \n6.5 | \n23.9 | \n423 | \n
| Koh Kong | \n5.4 | \n3.8 | \n0.4 | \n3.8 | \n6.2 | \n41 | \n
| Kratie | \n13.7 | \n9.1 | \n2.8 | \n10.2 | \n15.9 | \n138 | \n
| Mondul Kiri | \n42.7 | \n19.0 | \n3.6 | \n19.4 | \n45.2 | \n34 | \n
| Phnom Penh | \n10.3 | \n5.1 | \n1.3 | \n5.8 | \n11.0 | \n904 | \n
| Preah Vihear | \n38.1 | \n20.2 | \n8.5 | \n22.4 | \n40.0 | \n108 | \n
| Prey Veng | \n25.0 | \n7.2 | \n1.7 | \n8.6 | \n26.3 | \n366 | \n
| Pursat | \n10.8 | \n3.6 | \n0.7 | \n3.6 | \n10.8 | \n111 | \n
| Ratanak Kiri | \n24.7 | \n6.6 | \n0.8 | \n7.1 | \n25.7 | \n96 | \n
| Siemreap | \n21.5 | \n12.0 | \n2.9 | \n13.5 | \n24.6 | \n482 | \n
| Preah Sihanouk | \n10.3 | \n3.4 | \n1.6 | \n4.0 | \n10.8 | \n75 | \n
| Stung Treng | \n29.3 | \n12.5 | \n7.5 | \n16.3 | \n31.5 | \n66 | \n
| Svay Rieng | \n10.9 | \n8.3 | \n2.1 | \n8.8 | \n13.6 | \n217 | \n
| Takeo | \n18.4 | \n7.4 | \n1.2 | \n8.0 | \n19.3 | \n344 | \n
| Otdar Meanchey | \n15.1 | \n6.1 | \n0.0 | \n6.1 | \n15.4 | \n77 | \n
| Kep | \n8.0 | \n7.2 | \n1.2 | \n7.7 | \n11.0 | \n18 | \n
| Pailin | \n11.3 | \n9.5 | \n1.7 | \n9.5 | \n14.2 | \n29 | \n
| Tboung Khmum | \n11.1 | \n7.3 | \n2.3 | \n7.3 | \n14.0 | \n252 | \n
Now, we load geographical data of the administrative division of Cambodia. These data are obtained from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA). The information comes as a shapefile .shp and is available here in a zip file named som_adm_ocha_20250108_AB_SHP.zip. I downloaded this file and located it in /khmadm_gov_20181004_ab_shp folder.
shp_cambodia <- read_sf(\"khmadm_gov_20181004_ab_shp/khm_admbnda_adm1_gov_20181004.shp\")Next step is to combine the geographic information from the shapefile, shp_cambodia, with the data frame with the information on IPV, ipv_cambodia. To achieve this goal, we need to make sure that both datasets have a common variable to merge them. In this case, we will use the province name ADM1_EN in the shapefile and Region in the data frame. Then, I compare if the names of the provinces in both datasets match:
setdiff(shp_cambodia$ADM1_EN, ipv_cambodia$Region)[1] \"Oddar Meanchey\"The output indicates that there is one discrepancy. In shp_cambodia, the name of the province \"Oddar Meanchey\" corresponds to \"Otdar Meanchey\" in ipv_cambodia. I will rename this province in the shapefile to match the name in the data from the official report.
After this modification, I join both datasets using a left join to keep all the provinces in the shapefile and add the corresponding IPV values from the data frame:
\n\nMy idea is to develop three maps, each representing one type of violence against women and girls: emotional, physical, and sexual. I will first create the map corresponding to emotional violence.
\n# More font alternatives can be found here https://fonts.google.com/\nfont_add_google(\"Scada\", \"Scada\") \n\nshowtext_auto()\n\nmap_ipv_emotional_cambodia <- ggplot() +\n geom_sf(\n data = shp_cambodia, \n aes(fill = `Emotional violence`), \n color = \"white\", # Color of the province borders\n linewidth = 0.5\n ) +\n geom_sf_text(\n data = shp_cambodia,\n aes(label = `Emotional violence`),\n check_overlap = TRUE, # To avoid overlaping text in the plot\n size = 4, \n color = \"black\",\n fontface = \"bold\"\n ) +\n scale_fill_gradient2(\n low = \"#FFFE9EB3\",\n mid = \"#C53270B3\",\n high = \"#040404B3\",\n midpoint = 22.5,\n limits = c(0, 45),\n breaks = c(0, 15, 30, 45),\n labels = c(\"0\", \"15\", \"30\", \"45\")\n ) +\n labs(\n title = \"Emotional\", # Title of this subplot\n subtitle = \"\",\n caption = \"\",\n x = \"\",\n y = \"\",\n fill = \"Percentage of victims (%)\" # Title of the legend\n ) +\n theme_minimal(\n base_family = \"Scada\"\n ) +\n theme(\n plot.margin = margin(0, 40, 0, 40),\n \n plot.background = element_blank(), # No color in the background\n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 22,\n hjust = 0.5,\n margin = margin(5, 0, 5, 0)\n )\n )Now, I replicate the same format for physical and sexual violence.
\nmap_ipv_physical_cambodia <- ggplot() +\n geom_sf(\n data = shp_cambodia, \n aes(fill = `Physical violence`), \n color = \"white\", \n linewidth = 0.5\n ) +\n geom_sf_text(\n data = shp_cambodia,\n aes(label = `Physical violence`),\n check_overlap = TRUE,\n size = 4, \n color = \"black\",\n fontface = \"bold\"\n ) +\n scale_fill_gradient2(\n low = \"#FFFE9EB3\",\n mid = \"#C53270B3\",\n high = \"#040404B3\",\n midpoint = 22.5,\n limits = c(0, 45)\n ) +\n labs(\n title = \"Physical\",\n subtitle = \"\",\n caption = \"\",\n x = \"\",\n y = \"\",\n fill = \"\"\n ) + \n guides(fill = \"none\") +\n theme_minimal(\n base_family = \"Scada\"\n ) +\n theme(\n plot.margin = margin(0, 40, 0, 40),\n \n plot.background = element_blank(),\n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 22,\n hjust = 0.5,\n margin = margin(5, 0, 5, 0)\n )\n )\n\nmap_ipv_sexual_cambodia <- ggplot() +\n geom_sf(\n data = shp_cambodia, \n aes(fill = `Sexual violence`), \n color = \"white\", \n linewidth = 0.5\n ) +\n geom_sf_text(\n data = shp_cambodia,\n aes(label = `Sexual violence`),\n check_overlap = TRUE,\n size = 4, \n color = \"black\",\n fontface = \"bold\"\n ) +\n scale_fill_gradient2(\n low = \"#FFFE9EB3\",\n mid = \"#C53270B3\",\n high = \"#040404B3\",\n midpoint = 22.5,\n limits = c(0, 45)\n ) +\n labs(\n title = \"Sexual\",\n subtitle = \"\",\n caption = \"\",\n x = \"\",\n y = \"\",\n fill = \"\"\n ) + \n guides(fill = \"none\") +\n theme_minimal(\n base_family = \"Scada\"\n ) +\n theme(\n plot.margin = margin(0, 40, 0, 40),\n \n plot.background = element_blank(),\n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 22,\n hjust = 0.5,\n margin = margin(5, 0, 5, 0)\n )\n )Now, using {patchwork} package, I combine the three plots into one single representation. I also use plot_annotation() function to add a title, subtitle, and caption to the overall design.
# Design display two rows, one with the three maps, and the second one only with the common guide\nmap_ipv_cambodia <- (map_ipv_emotional_cambodia + map_ipv_physical_cambodia + map_ipv_sexual_cambodia) / guide_area() + \n # Custom title, subtitle, and caption\n plot_annotation(\n title = \"Intimate partner violence in Cambodia, 2021-22\",\n subtitle = \"Proportion of ever-partnered women and girls aged 15-49 who have been victims of violence perpetrated by their intimate partner, by type of violence.\",\n caption = \"**Data:** 2021-22 Cambodia Demographic and Health Survey (CDHS) <br> **Chart:** Juan Torres Mungu\u00eda\",\n theme = theme(\n # Title settings\n plot.title.position = \"plot\", \n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 30,\n hjust = 0,\n margin = margin(5, 0, 5, 50)\n ),\n \n plot.subtitle = element_text(\n color = \"grey20\",\n size = 25,\n hjust = 0,\n margin = margin(5, 0, 40, 50)\n ),\n \n # Caption settings\n # I use element_markdown() because the caption has markdown text\n plot.caption = element_markdown(\n color = \"black\",\n size = 15,\n hjust = 0,\n margin = margin(80, 0, 5, 50)\n )\n )) +\n\n plot_layout(guides = \"collect\",\n heights = unit(c(13, 1),\n c('cm', 'cm'))\n ) + \n theme_minimal(\n base_family = \"Scada\"\n ) &\n # Custom theme settings\n theme(\n legend.background = element_blank(),\n legend.position = \"bottom\",\n legend.title = element_text(\n margin = margin(5, 0, 10, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 14,\n hjust = 0.5\n ),\n legend.title.position = \"top\",\n legend.text = element_text(\n margin = margin(5, 5, 5, 5), \n face = \"bold\",\n color = \"grey10\",\n size = 14\n ),\n legend.key.width = unit(1.4, \"cm\"),\n legend.key.height = unit(1, \"cm\"),\n legend.direction = \"horizontal\",\n\n plot.background = element_rect(\n color = \"#E5E5E5\",\n fill = \"#E5E5E5\"),\n # Axis settings\n axis.title = element_blank(),\n axis.line = element_blank(),\n axis.text = element_blank(),\n axis.ticks = element_blank(),\n panel.grid = element_blank()\n \n )Finally, I use grid.lines() from {grid} package to add a horizontal line after the subtitle.
map_ipv_cambodia\n\ngrid.lines(\n x = c(0.04, 0.96), # Values should be between 0 and 1\n y = 0.88, # Value should be between 0 and 1\n gp = gpar(col = \"#C53270B3\", \n lwd = 4)\n)
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Three Maps to Examine the Prevalence of Intimate Partner\n Violence {(IPV)} in {Cambodia}},\n date = {2025-05-10},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/map-intimate-partner-violence-cambodia/},\n langid = {en}\n}\nUnlike some other agencies we track, NOAA grants lack an agency-specific data source for grants.\nFor NOAA, we rely primarily on USAspending, DOGE, and crowd-sourced information to make determinations on grant status.
\nDOGE provides data on terminated awards, although we have found that it is incomplete or occasionally incorrect, so we generally do manual review of DOGE data or use it only in combination with other sources.
\nWe use this source primarily for metadata including grantee information, grant assistance listing (CFDA), grant notice of funding announcement (NOFO), and total obligated and total outlaid amounts.\nWe also track changes to project end dates (i.e. project end date moved earlier, a \"cutoff\", or later, a \"reextension\"), and deobligation events (i.e. a negative amount obligated) which we use in combination as signals of disruption.
\nFile C is released roughly once a month (except that Sep/Oct are combined into one release).\nWe use this source primarily for outlay actions which we use to confirm reinstatements.
\nAn award is marked as terminated if it has been reported as such directly to us by a PI or via a trusted source, if it was listed by DOGE as being terminated, or if we see both a project end date cutoff and a deobligation in USAspending data.
\nAn award is labeled as reinstated if it has been reported as such directly to us by a PI, or if there is a re-extension of project end date or a positive obligation after a deobligation in USAspending awards data.\nReinstatements are confirmed by verifying that an award has received outlays since the reinstatement.
\nThe funds obligated and spent at time of earliest disruption are calculated using the current total obligated amount minus obligations since the earliest termination date and the total outlays minus outlays since the earliest termination date, respectively.\nThe outlay and obligation events come from File C, while the current total obligated and outlaid numbers are from the more regularly updating USAspending awards data.\nAs a result, these estimates may fluctuate week-to-week.
","doi":"https://doi.org/10.59350/74brw-em383","guid":"https://doi.org/10.59350/74brw-em383","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1790553600,"rid":"c2qz2-pwj80","summary":"Unlike some other agencies we track, NOAA grants lack an agency-specific data source for grants. For NOAA, we rely primarily on USAspending, DOGE, and crowd-sourced information to make determinations on grant status.","tags":["Methods","National Oceanic And Atmospheric Administration"],"title":"NOAA Grant Disruption Methodology","updated_at":1790621962,"url":"https://grantwitness.org/noaa/updates/methodology","version":"v1"},{"authors":[{"contributor_roles":[],"family":"Shan","given":"Mridula","url":"https://orcid.org/0000-0003-4613-2669"},{"affiliation":[{"id":"https://ror.org/019jywm96","name":"rOpenSci"}],"contributor_roles":[],"family":"Ross","given":"Noam","url":"https://orcid.org/0000-0002-2136-0000"}],"blog":{"authors":null,"community_id":"3446712e-ab03-4bb7-8641-571dcda8d8cd","created":1787702400,"current_feed_url":null,"description":"Tracking disruptions to U.S. federal science funding, grant by grant.","doi":"https://doi.org/10.59350/grantwitness","favicon":"https://rogue-scholar.org/api/communities/3446712e-ab03-4bb7-8641-571dcda8d8cd/logo","feed_format":"application/atom+xml","feed_url":"https://grantwitness.org/updates/rogue-scholar.xml","filter":null,"generator":"Other","home_page_url":"https://grantwitness.org/updates?type=updates","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"grantwitness","status":"active","subfield":"3321","title":"Grant Witness","updated":1790553600,"use_api":null},"blog_name":"Grant Witness","blog_slug":"grantwitness","content_html":"At Grant Witness, we have been tracking the damage of funding terminations across federal agencies since March 2025. We have focused on our core research agencies (NIH and NSF) and expanded to other agencies in science and public health such as EPA and CDC. Yet scientific research, development, and application are core functions across the federal government, and we are always striving to expand our view to better understand how the current administration is upending the American science enterprise.
\nToday we're expanding both our termination and grantmaking trackers to include the National Oceanic and Atmospheric Administration (NOAA). NOAA is a scientific agency within the United States Department of Commerce. Its stated mission is \"to understand and predict our changing environment, from the deep sea to outer space, and to manage and conserve America's coastal marine resources.\" It includes the National Weather Service.
\nNOAA's work has been targeted from the start of the administration. The Trump administration has engaged in staff firings, suppression of scientific reports, and proposed budget cuts of more than 25% for both FY 2026 and 2027.
\nLike many other agencies, NOAA terminated research and support grants en masse starting in April 2025 under the disastrous tenure of the \"Department of Government Efficiency\" (DOGE). Neither NOAA nor DOGE has published a reliable accounting of these terminations. Grant Witness, through reporting and data mining, has identified 26 NOAA grants that have been terminated, only 4 of which have been reinstated.
\nThese grants span 13 government programs (CFDAs, Table\u00a01) across 4 of the 6 line offices at NOAA (Table\u00a02). Together, they represent a loss of $28,237,466 in funding to states, cities, nonprofits and universities. Climate research is the most affected area, but a variety of programs in coastal and marine management have been shuttered.
\n| CFDA Program Title | \nDisrupted Grants | \nAmt Remaining at Disruption | \n
|---|---|---|
| Climate and Atmospheric Research | \n11 | \n$11,217,919 | \n
| Cooperative Institute Inter Agency Funded Activities | \n2 | \n$1,927,901 | \n
| Marine Sanctuary Program | \n2 | \n$235,000 | \n
| Office for Coastal Management | \n2 | \n$1,470,012 | \n
| Chesapeake Bay Studies | \n1 | \n$40,700 | \n
| Coastal Zone Management Administration Awards | \n1 | \n$179,600 | \n
| Coastal Zone Management Estuarine Research Reserves | \n1 | \n$8,486.71 | \n
| Cooperative Science and Education Program | \n1 | \n$11,600 | \n
| Habitat Conservation | \n1 | \n$15.51 | \n
| National Oceanic and Atmospheric Administration Noaa Cooperative Institutes | \n1 | \n$2,757,930 | \n
| Noaa Mission Related Education Awards | \n1 | \n$364,360 | \n
| Pacific Fisheries Data Program | \n1 | \n$766,710 | \n
| Sea Grant Support | \n1 | \n$9,257,231 | \n
| Line Office Name | \nDisrupted Grants | \nAmt Remaining at Disruption | \n
|---|---|---|
| Office of Oceanic and Atmospheric Research | \n15 | \n$25,160,981 | \n
| National Ocean Service | \n6 | \n$1,893,099 | \n
| National Marine Fisheries Service | \n4 | \n$819,026 | \n
| Office of the Under Secretary | \n1 | \n$364,360 | \n
As one would expect, coastal states have borne the brunt of these terminations. New York has the greatest number of disrupted NOAA grants followed by Washington, Massachusetts, and New Jersey (tied at three disrupted grants) (Table\u00a03). The states with the greatest loss in funding include Washington, North Crolina, and New Jersey.
\n| State | \nDisrupted Grants | \nRemaining at Disruption | \n
|---|---|---|
| New York | \n4 | \n$787,780 | \n
| Washington | \n3 | \n$11,918,505 | \n
| New Jersey | \n3 | \n$4,685,831 | \n
| Massachusetts | \n3 | \n$384,447 | \n
| North Carolina | \n2 | \n$5,725,610 | \n
| Virginia | \n2 | \n$124,667 | \n
| California | \n1 | \n$1,921,064 | \n
| Louisiana | \n1 | \n$1,126,012 | \n
| Oregon | \n1 | \n$766,710 | \n
| Rhode Island | \n1 | \n$344,000 | \n
| Wyoming | \n1 | \n$271,297 | \n
| Michigan | \n1 | \n$89,519 | \n
| Illinois | \n1 | \n$75,000 | \n
| Alaska | \n1 | \n$17,008 | \n
| Maine | \n1 | \n$16 | \n
On our NOAA agency pulse page, we are tracking the new awards made by the agency. Unlike many of the agencies we track, NOAA does not provide up-to-date grant listings, so we use USAspending data to track new awards, which can lag by over a month. Nonetheless, we can see NOAA is far behind in its distribution of Congressionally-appropriated funds. Unlike FY 2025, when many NOAA grants were distributed in the closing days of the Biden administration, NOAA was slow to start this year and has not made new grants since mid-August. It has distributed 61% fewer funds and 74% fewer grants than the 2021\u20132024 average. We will see if the agency can catch up in the remaining weeks of the fiscal year.
\n
We will continue to track the changes and disruptions at NOAA. Keep an eye on the agency page, or follow our updates or newsletter. As with all our work at Grant Witness, crowdsourced information is key to building and verifying our databases. Please continue to submit reports of grant terminations or reinstatements using our web form, and don't hesitate to contact us.
","doi":"https://doi.org/10.59350/zbrpd-kbv50","guid":"https://doi.org/10.59350/zbrpd-kbv50","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1790553600,"rid":"cyhts-3s181","summary":"At Grant Witness, we have been tracking the damage of funding terminations across federal agencies since March 2025. We have focused on our core research agencies (NIH and NSF) and expanded to other agencies in science and public health such as EPA and CDC.","tags":["Update","National Oceanic And Atmospheric Administration"],"title":"Grant Witness is tracking NOAA Terminations and Funding","updated_at":1790621961,"url":"https://grantwitness.org/noaa/updates/2026-08-24-NOAA-targeted-ncrs","version":"v1"},{"authors":[{"contributor_roles":[],"family":"Jivkova","given":"Katrin"}],"blog":{"authors":null,"community_id":"1e54953c-587a-4743-89af-ebd2cb45660c","created":1789516800,"current_feed_url":null,"description":null,"doi":"https://doi.org/10.59350/pensoft","favicon":"https://rogue-scholar.org/api/communities/1e54953c-587a-4743-89af-ebd2cb45660c/logo","feed_format":"application/atom+xml","feed_url":"https://blog.pensoft.net/feed/atom/","filter":null,"generator":"WordPress","home_page_url":"https://blog.pensoft.net/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"pensoft","status":"active","subfield":"2303","title":"Pensoft Blog","updated":1790604504,"use_api":true},"blog_name":"Pensoft Blog","blog_slug":"pensoft","content_html":"\n\nPensoft Publishers is pleased to announce the publication of the first issue of the Journal of Scientific and Technological Innovations (STI), an international journal for cutting-edge research across the scientific and technological spectrum.
\n\n\n\nThe journal notably adopts Pensoft's ARPHA platform, an innovative end-to-end full-featured publishing solution that streamlines the entire publication workflow \u2013 from initial manuscript submission to indexing, distribution, hosting and archiving through one coordinated platform. To date, ARPHA provides \"a home\" to about 100 open-access journals across scientific fields.
\n\n\nThe Journal of Scientific and Technological Innovations is an international, multidisciplinary academic journal supervised by Ho Chi Minh City University of Industry and Trade (HUIT). The journal publishes original, peer-reviewed research contributions.
\n\n\n\nIts editorial purview encompasses a broad range of scientific and technological domains:
\n\n\n\nThe journal is led by an international team of scholars:
\n\n\n\nEditor-in-Chief: Thi Hong Anh Le, Ho Chi Minh City University of Industry and Trade, Vietnam
\n\n\n\nDeputy Editors-in-Chief: Tuyen Truong, RMIT University, Australia; Thanh-Binh Nguyen, National Kaohsiung University of Science and Technology, Taiwan
\n\n\n\nThe first issue of STI presents six contributions spanning food science, materials science, and applied food technology:
\n\n\n\nA cost-benefit analysis of small-holder floating-house fish farming in Chau Doc City, Vietnamese Mekong Delta
\n\n\n\nA study of 35 Mekong Delta fish farmers reveals significant disparities in farm profitability across different river reaches, with operations facing demographic vulnerabilities and heavy reliance on unpaid labour.
\n\n\n\nRead more: Farzana H, Tran DD, Nguyen TD, Park E (2026) A cost-benefit analysis of small-holder floating-house fish farming in Chau Doc City, Vietnamese Mekong Delta. Journal of Scientific and Technological Innovations 1: e209236. https://doi.org/10.3897/sti.e209236
\n\n\n\nUltrasound-assisted pullulanase modification increases resistant starch content and alters the structure\u2013function properties of water caltrop starch
\n\n\n\nA combined ultrasound and enzymatic treatment increased resistant starch content in water caltrop starch from 21.72% to 28.15%, whilst enhancing its water and oil absorption capacities. The approach offers potential for developing resistant starch-enriched functional food ingredients.
\n\n\n\nRead more: Phan VM, Ngo TMP, Do TTN, Luyen TTH, Duong TT, Tran NM, Tran CH (2026) Ultrasound-assisted pullulanase modification increases resistant starch content and alters the structure\u2013function properties of water caltrop starch. Journal of Scientific and Technological Innovations 1: e209539. https://doi.org/10.3897/sti.e209539
\n\n\n\nEffects of Melatonin/\u03b2-cyclodextrin Treatment on Physicochemical and Bioactive Attributes of Cherry Tomatoes during Storage
\n\n\n\nMelatonin complexed with \u03b2-cyclodextrin reduced cherry tomato weight loss by up to 56% during refrigerated storage and proved more effective than free melatonin in preserving colour, acidity, and polyphenol content. The complex formation provides a viable strategy for extending shelf-life in harvested produce.
\n\n\n\nRead more: Tran THC, Dong TT, Nguyen PH,Le NH,Kuo C-H (2026) Effects of Melatonin/\u03b2-cyclodextrin Treatment on Physicochemical and Bioactive Attributes of Cherry Tomatoes during Storage. Journal of Scientific and Technological Innovations 1: e212043. https://doi.org/10.3897/sti.e212043
\n\n\n\nDesign, Functionalities, and Emerging Applications of Ti3C2Tx MXene-Based Composites: A Review
\n\n\n\nThis review synthesises current knowledge on titanium carbide-based MXenes, two-dimensional materials with exceptional conductivity and mechanical properties. MXene-based composites demonstrate applications in gas sensing, photocatalytic carbon dioxide reduction, and energy storage devices.
\n\n\n\nRead more: Ta QTH, Pham NN, Tran QN, Lu TMT, Noh JS (2026) Design, Functionalities, and Emerging Applications of Ti3C2Tx MXene-Based Composites: A Review. Journal of Scientific and Technological Innovations 1: e205045. https://doi.org/10.3897/sti.\u0435205045
\n\n\n\nChemical Composition, Bioactivities, and Cosmetic Applications of Coconut Oil and Its Derivatives
\n\n\n\nRecent evidence demonstrates that coconut oil, when appropriately applied, offers dermatological benefits, with antimicrobial, anti-inflammatory, and antioxidant properties. The research clarifies interactions between coconut oil and skin microflora and explains previous inconsistencies in published findings.
\n\n\n\nRead more: Dat NM, Phong DT, Long HB, Dat NT, Thy LTM (2026) Chemical Composition, Bioactivities, and Cosmetic Applications of Coconut Oil and Its Derivatives. Journal of Scientific and Technological Innovations 1: e205088. https://doi.org/10.3897/sti.e205088
\n\n\n\nRetention of bilberry anthocyanins in wood hemicelluloses coated microcapsules under the impact of humidity and packaging conditions
\n\n\n\nStorage at 8% relative humidity limited bilberry anthocyanin degradation to 17-22% over 24 weeks, with wood hemicelluloses providing comparable protection to conventional gum Arabic. The findings establish moisture control as important for developing shelf-stable natural functional powders with sustainable coating materials.
\n\n\n\nThe diversity of contributions in this first issue reflects STI’s multidisciplinary scope. Researchers interested in submitting manuscripts are invited to consult the journal's guidelines and submission portal.
\n\n\n\n***
\n\n\n\nAccess the full first issue here.
\n\n","doi":"https://doi.org/10.59350/12ke2-9s798","guid":"https://blog.pensoft.net/?p=21014","image":"https://blog.pensoft.net/wp-content/uploads/2026/09/SciComm-PR-2026-Instagram-2026-09-28T141749.018-1.png","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1790553600,"rid":"t1tqd-psc54","summary":"Tweet Pensoft Publishers is pleased to announce the publication of the first issue of the Journal of Scientific and Technological Innovations (STI), an international journal for cutting-edge research across the scientific and technological spectrum.","tags":["Academic Publishing","ARPHA Platform","Interdisciplinary Research","New Journal","Scientific Innovation"],"title":"New Journal of Scientific and Technological Innovations launches on ARPHA","updated_at":1790604910,"url":"https://blog.pensoft.net/2026/09/28/new-journal-of-scientific-and-technological-innovations-launches-on-arpha/","version":"v1"},{"authors":[{"contributor_roles":[],"family":"Torres Munguia","given":"Juan Armando","url":"https://orcid.org/0000-0003-3432-6941"}],"blog":{"authors":null,"community_id":"ae893dd4-41e3-426b-8b48-6c9c75639805","created":1761523200,"current_feed_url":null,"description":"Exploring the art and science of working with development and humanitarian data","doi":"https://doi.org/10.59350/juan_torresmunguia","favicon":"https://rogue-scholar.org/api/communities/ae893dd4-41e3-426b-8b48-6c9c75639805/logo","feed_format":"application/rss+xml","feed_url":"https://juan-torresmunguia.netlify.app/blog/index.xml","filter":null,"generator":"Quarto","home_page_url":"https://juan-torresmunguia.netlify.app/blog/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"juan_torresmunguia","status":"active","subfield":"1804","title":"Data Art & Science Blog","updated":1783172007,"use_api":null},"blog_name":"Data Art & Science Blog","blog_slug":"juan_torresmunguia","content_html":"In these notes, I will show you how to use {ggplot} to create a waffle chart about the frequency of disease outbreaks in the world. The final output will look like this:
\n
The source of information is the global dataset of pandemic- and epidemic-prone disease outbreaks, whose data is freely available at the GitHub repository of the disease outbreaks project.
\nThe unit of analysis in the database is an outbreak, defined as the occurrence of at least one case of a specific disease in a country -or territory- during a particular year. Therefore, a country -or territory- cannot have more than one outbreak of the same disease in the same year, although it may experience outbreaks of different diseases within the same year. A country can only have multiple outbreaks of the same disease if they occur in different years.
\nThe dataset includes information on over 3,000 pandemic- and epidemic-prone disease outbreaks associated with more than 80 different infectious diseases, occurring globally from January 1996 to the present. The diseases are classified according to the International Classification of Diseases, 10th Revision (ICD-10), and the dataset contains information on the year, country, and disease for each outbreak.
\nThe dataset of pandemic- and epidemic-prone disease outbreaks is also part of the Humanitarian Data Exchange coordinated by the United Nations Office for the Coordination of Humanitarian Affairs (OCHA).
\nTo create the waffle chart, we will use the following R packages:
\n\nThe dataset is stored in an Excel. We will download it from the GitHub repository and load it into R:
\n# Define the GitHub raw content URL\nlink_data <- \"https://github.com/jatorresmunguia/disease_outbreak_news/raw/refs/heads/main/Humanitarian%20Data%20Exchange/data_2_share/disease_outbreaks_HDX.xlsx\"\n\n# Create a temporary file to store the data\ntemp <- tempfile(fileext = \".xlsx\")\n\n# Download the data from GitHub using authentication\nget_req <- GET(\n link_data,\n authenticate(Sys.getenv(\"GITHUB_PAT\"), \"\"),\n write_disk(path = temp)\n)\n\n# Read the Excel file into R\noutbreaks <- read_excel(temp)We will group the data by year and disease, and then filter the top 10 diseases from 1996 to 2024. We will also group the data by year and disease, and then filter the top 10 diseases from 2015 to 2024. Finally, we will create a factor variable for the top 10 diseases and the \"Other\" category:
\n# Top 10 diseases 1996-2024\ntop_10_icd104n <- outbreaks |>\n filter(Year < 2025) |>\n group_by(icd104n) |>\n summarise(freq = n(), .groups = \"drop\") |>\n arrange(-freq) |>\n slice_head(n = 10)\n\n# Only top-10 categories and Others\noutbreaks_year_disease_grouped <- outbreaks |>\n filter(Year < 2025 & Year > 2014) |>\n mutate(icd104n = ifelse(icd104n %in% top_10_icd104n$icd104n, icd104n, \"Other\")) |> # Group the diseases\n group_by(Year, icd104n) |>\n summarise(freq = n(), .groups = \"drop\") |> # Count the number of outbreaks\n arrange(-freq)\n\n# Create a factor variable for the top 10 diseases and the \"Other\" category\noutbreaks_year_disease_grouped <- outbreaks_year_disease_grouped |>\n mutate(icd104n = factor(icd104n,\n levels = c(\n \"Acute poliomyelitis, unspecified\", \"Classical cholera\", \"COVID-19, virus identified\",\n \"Dengue, unspecified\", \"Influenza due to identified zoonotic or pandemic influenza virus\", \"Measles\",\n \"Meningococcal meningitis\", \"Middle East respiratory syndrome coronavirus [MERS-CoV]\", \"Monkeypox\",\n \"Yellow fever, unspecified\", \"Other\"\n ),\n labels = c(\n \"Acute poliomyelitis\", \"Cholera\", \"COVID-19\", \"Dengue\", \"Influenza\",\n \"Measles\", \"Meningococcal meningitis\", \"MERS-CoV\", \"Monkeypox\",\n \"Yellow fever\", \"Other\"\n )\n ))# Add custom font\nfont_add_google(\"Roboto Condensed\", \"Roboto Condensed\") \n\nshowtext_auto()\n\n# Custom theme for the waffle chart\ntheme_waffle_chart <- function() {\n theme_minimal(\n base_family = \"Roboto Condensed\" # Base theme with custom font\n ) +\n # Custom theme settings\n theme(\n # Axis settings\n axis.title = element_blank(), # Remove axis title\n axis.line = element_blank(), # Remove axis line\n \n # Title settings\n plot.title.position = \"plot\", # Position of the title\n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 24,\n margin = margin(5, 0, 5, 0), # Top, right, bottom, left\n width = unit(1, \"npc\") # Width of the title, npc == 1 corresponds to the full width of the plot\n ),\n \n # Subtitle settings\n plot.subtitle = element_textbox(\n color = \"grey50\",\n face = \"bold\",\n size = 14,\n margin = margin(0, 0, 10, 0),\n width = unit(1, \"npc\")\n ),\n \n # Legend settings\n legend.position = \"top\",\n legend.title = element_blank(),\n legend.spacing.x = unit(0.2, \"cm\"),\n legend.key.spacing = unit(0.5, \"cm\"), # Spacing between legend keys\n legend.text = element_text(\n margin = margin(5, 2, 5, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 9\n ),\n legend.direction = \"horizontal\",\n legend.byrow = FALSE,\n \n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n size = 8,\n hjust = 0\n ),\n plot.caption.position = \"plot\",\n\n plot.background = element_rect(\n color = \"white\",\n fill = \"white\"\n ),\n plot.margin = margin(20, 40, 20, 40)\n )\n}\n\n# Title, subtitle, and caption for the waffle chart\ntitle_chart <- \"Pandemic- and epidemic-prone disease outbreaks in the world\"\nsubtitle_chart <- \"In 2024, fewer countries reported public health events linked to COVID-19, while outbreaks of viral diseases transmitted by vectors increased\"\ncaption_chart <- paste0(\"**Note:** Each square represents a country.\",\n \"<br>\", \n \"**Data:** A global dataset of pandemic- and epidemic-prone disease outbreaks (DOI: 10.1038/s41597-022-01797-2)\",\n \"<br>\", \n \"**Graphic:** Juan Torres Mungu\u00eda\"\n )\n\n# Define the colors to be used in the waffle chart\ncolors <- c(\n \"#B66DFF\", \"#7e0021\", \"#ff420e\", \"#579d1c\", \"#83caff\",\n \"#0084d1\", \"#aecf00\", \"#4b1f6f\", \"#FF86FF\", \"#ffd320\", \"grey80\"\n )ggplot(\n outbreaks_year_disease_grouped,\n aes(fill = icd104n, values = freq)\n) +\n geom_waffle(\n size = 0.25, # Size of the waffle squares\n n_rows = 10, # Number of rows in the waffle chart\n flip = TRUE, # Flip the waffle chart\n color = \"white\", # Color of the border of the waffle squares\n make_proportional = FALSE\n ) + # Do not make the waffle chart proportional\n facet_wrap(~Year, # Facet the waffle chart by year\n nrow = 1,\n strip.position = \"bottom\"\n ) + # Position of the facet strip\n scale_fill_manual(values = colors) +\n scale_x_discrete() +\n scale_y_continuous(\n labels = function(x) x * 10,\n expand = c(0, 0)\n ) + # Expand the y-axis\n coord_equal() +\n labs(\n title = title_chart,\n subtitle = subtitle_chart,\n caption = caption_chart,\n x = \"\",\n y = \"Total number of outbreaks\",\n fill = \"\"\n ) +\n guides(fill = guide_legend(nrow = 2 # Number of rows in the legend\n )) + \n theme_waffle_chart()showtext_opts(dpi = 320) # Set the resolution of the image 320 dpi is for high-quality images (\"retina\")\nggsave(\n \"waffle-pandemics.png\",\n dpi = 320,\n width = 14,\n height = 7,\n units = \"in\"\n)\nshowtext_auto(FALSE)@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {A Waffle Chart of Disease Outbreaks in the World},\n date = {2025-01-12},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/waffle-chart-disease-outbreaks/},\n langid = {en}\n}\nIn this article, I will guide you through the process of analyzing Mexico's 2024 National Survey of Victimization and Perception of Public Safety (Encuesta Nacional de Victimizaci\u00f3n y Percepci\u00f3n sobre Seguridad P\u00fablica, ENVIPE) using the {tidyverse} framework for survey analysis with the {srvyr} package. Specifically, we will calculate the crime prevalence rate per 100,000 inhabitants by state, according to the victim's sex ( Tasa de prevalencia delictiva por entidad federativa por cada cien mil habitantes, seg\u00fan sexo de la v\u00edctima).
\nOfficial tables with this indicator can be found here.
\nBy the end of this article, you will be able to use the {tidyverse} to:
\nENVIPE is an official survey that provides nationally and state-representative information on household crime victimization and perceptions of public safety in Mexico. Conducted annually since 2011 by the National Institute of Statistics and Geography (INEGI), ENVIPE collects data on the following crimes:
\nAll survey-related information, including questionnaires, methodology, glossaries, summaries, and microdata, is fully available at INEGI's website.
\nFirst, install and load the required packages:
\n\nMicrodata is available here in various formats (CSV, DBF, DTA, SAV, and RData). In this post, we will use the RData files.
\nload(archive_read(\n archive = \"https://www.inegi.org.mx/contenidos/programas/envipe/2024/microdatos/bd_envipe_2024_RData.zip\",\n file = \"BD_ENVIPE_2024.RData\"\n))Using the {archive} package, we can directly extract the required file from bd_envipe_2024_RData.zip.
This loads the following objects into the global environment:
\nTHogar: This object contains the interview result, the row of the selected respondent, the household factor, and other household characteristics.
TMod_Vic: This object includes information on the crimes suffered by the selected respondent and their household during the reference year, specifically during 2023. Information comes from the victimization module of the survey.
TPer_Vic1: This object includes information related to the perception of public security within the respondent's geographic area, antisocial behaviors in their immediate surroundings, and changes in habits due to fear of becoming a crime victim in 2023. Additionally, it presents information regarding the performance of authorities and trust in them. In other words, this table covers Sections IV and V of the main questionnaire.
TPer_Vic2: This object consolidates information related to household and personal victimization, particularly focusing on high-impact crimes experienced by the household (such as kidnapping, enforced disappearance, or homicide) or by the selected respondent. This table includes Sections VI and VII of the main questionnaire.
TSDem: This object contains the socio-demographic characteristics of household members, corresponding to Section III of the main questionnaire.
TSVivienda: This object contains the characteristics of the households captured through the questionnaire cover page, as well as the count of household residents and the number of households. In other words, it includes Sections I and II of the main questionnaire.
To estimate the crime prevalence rate, we use information from TPer_Vic2 (screening questions on victimization from the main questionnaire) and TMod_Vic (specific victimization module applied to the respondents who reported being a victim of at least one of the crimes included in the survey during the reference period). TPer_Vic2 will serve to characterize the total population and TMod_Vic to describe the victims and the associated crimes.
ENVIPE follows a probabilistic, three-stage, stratified, and clustered sampling design. The sampling unit is the household, while the unit of analysis is the population aged 18 and older in selected households. More details on the sampling design can be found here(in Spanish).
\nFrom the documentation on data structure(in Spanish), the key variables for specifying the survey design are:
\nID_PER: Person identifier
EST_DIS: Stratum
UPM_DIS: Primary sampling unit
FAC_ELE: Weighting factor for public security perception and victimization estimates.
We specify the survey design as follows:
\n# Identify victims from TMod_Vic, excluding\n# \"Wall graffiti, intentional vehicle damage, and other vandalism\" (BPCOD != \"03\")\nid_victims <- TMod_Vic |>\n filter(BPCOD != \"03\") |>\n pull(ID_PER)\n\n# Create victimization variable in TPer_Vic2\nTPer_Vic2 <- TPer_Vic2 |>\n mutate(victimization = ifelse(ID_PER %in% id_victims, 1, 0))\n\n# Convert FAC_ELE (weights) to numeric\nTPer_Vic2 <- TPer_Vic2 |>\n mutate(FAC_ELE = as.numeric(FAC_ELE))\n\n# Label SEXO variable\nTPer_Vic2 <- TPer_Vic2 |>\n mutate(SEXO = recode(SEXO, \"1\" = \"Man\", \"2\" = \"Woman\"))\n\n# Create survey object\nenvipe_design <- TPer_Vic2 |>\n as_survey_design(\n weights = FAC_ELE, # Weights (inverse probability)\n strata = EST_DIS, # Stratum\n ids = UPM_DIS, # Primary sampling unit\n nest = TRUE # Nesting within strata\n )To compute proportions from survey data, we use survey_mean() from {srvyr} package. From ENVIPE's documentation, we know that NOM_ENT represents the state, while SEXO corresponds to the respondent's sex.
envipe_design |>\n # `interact` function is used to calculate the proportion\n # over the interaction of state and sex\n group_by(interact(NOM_ENT, SEXO)) |>\n summarize(\n Proportion = survey_mean(victimization,\n level = 0.90, # Confidence level\n vartype = c(\"se\", \"ci\", \"cv\") # Variability indicators\n )\n ) |>\n ungroup() |>\n mutate(\n `Prevalence rate` = Proportion * 100000, # Rate per 100,000 people\n `Standard error` = Proportion_se * 100000,\n `Low 90% CI` = Proportion_low * 100000,\n `Upper 90% CI` = Proportion_upp * 100000,\n `Coefficient of variation %` = Proportion_cv * 100,\n State = str_to_sentence(NOM_ENT) # Capitalize only the first letter\n ) |>\n rename(Sex = SEXO) |>\n select(\n State, Sex, `Prevalence rate`, `Standard error`,\n `Low 90% CI`, `Upper 90% CI`, `Coefficient of variation %`\n ) |>\n kbl() |> # HTML table\n kable_styling() # Apply bootstrap theme to the tableThe output table looks like this:
\n| State | \nSex | \nPrevalence rate | \nStandard error | \nLow 90% CI | \nUpper 90% CI | \nCoefficient of variation % | \n
|---|---|---|---|---|---|---|
| Aguascalientes | \nMan | \n33710.81 | \n2172.7372 | \n30136.71 | \n37284.91 | \n6.445224 | \n
| Aguascalientes | \nWoman | \n32041.57 | \n1831.7081 | \n29028.45 | \n35054.68 | \n5.716662 | \n
| Baja California | \nMan | \n23043.35 | \n1604.9769 | \n20403.20 | \n25683.50 | \n6.965032 | \n
| Baja California | \nWoman | \n24480.64 | \n1536.1771 | \n21953.66 | \n27007.61 | \n6.275070 | \n
| Baja California Sur | \nMan | \n19329.94 | \n1244.1353 | \n17283.36 | \n21376.51 | \n6.436314 | \n
| Baja California Sur | \nWoman | \n23817.58 | \n1466.2963 | \n21405.56 | \n26229.60 | \n6.156361 | \n
| Campeche | \nMan | \n21832.87 | \n1455.8508 | \n19438.03 | \n24227.71 | \n6.668161 | \n
| Campeche | \nWoman | \n23213.74 | \n1381.0933 | \n20941.88 | \n25485.61 | \n5.949465 | \n
| Chiapas | \nMan | \n16602.99 | \n986.0052 | \n14981.03 | \n18224.94 | \n5.938722 | \n
| Chiapas | \nWoman | \n12074.36 | \n978.7862 | \n10464.28 | \n13684.44 | \n8.106318 | \n
| Chihuahua | \nMan | \n22303.22 | \n1234.5065 | \n20272.49 | \n24333.96 | \n5.535103 | \n
| Chihuahua | \nWoman | \n22804.31 | \n1181.0327 | \n20861.54 | \n24747.08 | \n5.178990 | \n
| Ciudad De Mexico | \nMan | \n34562.66 | \n1089.2714 | \n32770.84 | \n36354.49 | \n3.151584 | \n
| Ciudad De Mexico | \nWoman | \n30746.76 | \n990.3528 | \n29117.66 | \n32375.87 | \n3.220999 | \n
| Coahuila De Zaragoza | \nMan | \n18179.70 | \n1249.1422 | \n16124.89 | \n20234.51 | \n6.871084 | \n
| Coahuila De Zaragoza | \nWoman | \n20447.19 | \n1238.7111 | \n18409.54 | \n22484.84 | \n6.058099 | \n
| Colima | \nMan | \n21106.85 | \n1593.5970 | \n18485.43 | \n23728.28 | \n7.550140 | \n
| Colima | \nWoman | \n21295.28 | \n1351.1455 | \n19072.68 | \n23517.88 | \n6.344813 | \n
| Durango | \nMan | \n15240.74 | \n1133.1868 | \n13376.67 | \n17104.80 | \n7.435250 | \n
| Durango | \nWoman | \n17401.18 | \n1166.0910 | \n15482.99 | \n19319.37 | \n6.701219 | \n
| Guanajuato | \nMan | \n20057.38 | \n1428.0981 | \n17708.20 | \n22406.57 | \n7.120063 | \n
| Guanajuato | \nWoman | \n18930.15 | \n1287.0811 | \n16812.93 | \n21047.37 | \n6.799107 | \n
| Guerrero | \nMan | \n16981.89 | \n1494.6554 | \n14523.22 | \n19440.56 | \n8.801466 | \n
| Guerrero | \nWoman | \n14572.44 | \n1238.6234 | \n12534.93 | \n16609.94 | \n8.499769 | \n
| Hidalgo | \nMan | \n22397.82 | \n1338.0474 | \n20196.77 | \n24598.88 | \n5.974006 | \n
| Hidalgo | \nWoman | \n20215.02 | \n1157.0294 | \n18311.73 | \n22118.30 | \n5.723613 | \n
| Jalisco | \nMan | \n24390.13 | \n1548.7916 | \n21842.40 | \n26937.85 | \n6.350076 | \n
| Jalisco | \nWoman | \n25141.90 | \n1438.5170 | \n22775.58 | \n27508.23 | \n5.721592 | \n
| Mexico | \nMan | \n32906.59 | \n1632.3183 | \n30221.47 | \n35591.71 | \n4.960460 | \n
| Mexico | \nWoman | \n33023.78 | \n1483.3471 | \n30583.71 | \n35463.85 | \n4.491755 | \n
| Michoacan De Ocampo | \nMan | \n14405.17 | \n1031.1125 | \n12709.01 | \n16101.32 | \n7.157934 | \n
| Michoacan De Ocampo | \nWoman | \n15485.20 | \n1036.3320 | \n13780.46 | \n17189.94 | \n6.692401 | \n
| Morelos | \nMan | \n26995.43 | \n1517.6054 | \n24499.00 | \n29491.85 | \n5.621713 | \n
| Morelos | \nWoman | \n24869.74 | \n1688.5969 | \n22092.04 | \n27647.44 | \n6.789765 | \n
| Nayarit | \nMan | \n16075.56 | \n1516.6195 | \n13580.76 | \n18570.36 | \n9.434318 | \n
| Nayarit | \nWoman | \n18639.61 | \n1597.1842 | \n16012.28 | \n21266.94 | \n8.568766 | \n
| Nuevo Leon | \nMan | \n23602.12 | \n1261.9498 | \n21526.25 | \n25678.00 | \n5.346764 | \n
| Nuevo Leon | \nWoman | \n21440.72 | \n1211.6372 | \n19447.60 | \n23433.83 | \n5.651104 | \n
| Oaxaca | \nMan | \n12940.56 | \n1158.2366 | \n11035.29 | \n14845.83 | \n8.950434 | \n
| Oaxaca | \nWoman | \n13548.61 | \n1017.7020 | \n11874.51 | \n15222.70 | \n7.511489 | \n
| Puebla | \nMan | \n26519.64 | \n1342.0473 | \n24312.00 | \n28727.27 | \n5.060579 | \n
| Puebla | \nWoman | \n24176.45 | \n1124.4656 | \n22326.73 | \n26026.17 | \n4.651078 | \n
| Queretaro | \nMan | \n30297.50 | \n1620.6534 | \n27631.57 | \n32963.44 | \n5.349132 | \n
| Queretaro | \nWoman | \n25746.31 | \n1260.3524 | \n23673.06 | \n27819.56 | \n4.895273 | \n
| Quintana Roo | \nMan | \n22592.04 | \n1384.4250 | \n20314.70 | \n24869.39 | \n6.127932 | \n
| Quintana Roo | \nWoman | \n21775.64 | \n1237.8095 | \n19739.48 | \n23811.81 | \n5.684376 | \n
| San Luis Potosi | \nMan | \n24019.02 | \n1874.7438 | \n20935.11 | \n27102.93 | \n7.805247 | \n
| San Luis Potosi | \nWoman | \n23227.97 | \n1491.3111 | \n20774.80 | \n25681.14 | \n6.420324 | \n
| Sinaloa | \nMan | \n21181.44 | \n1122.0160 | \n19335.76 | \n23027.13 | \n5.297165 | \n
| Sinaloa | \nWoman | \n20562.22 | \n1086.2790 | \n18775.32 | \n22349.12 | \n5.282887 | \n
| Sonora | \nMan | \n26864.66 | \n1816.2677 | \n23876.94 | \n29852.37 | \n6.760808 | \n
| Sonora | \nWoman | \n25742.16 | \n1902.0264 | \n22613.37 | \n28870.95 | \n7.388759 | \n
| Tabasco | \nMan | \n25881.72 | \n1354.6212 | \n23653.40 | \n28110.04 | \n5.233891 | \n
| Tabasco | \nWoman | \n25522.45 | \n1091.8488 | \n23726.38 | \n27318.51 | \n4.277994 | \n
| Tamaulipas | \nMan | \n17788.15 | \n1075.8038 | \n16018.47 | \n19557.82 | \n6.047869 | \n
| Tamaulipas | \nWoman | \n18368.67 | \n1017.1639 | \n16695.46 | \n20041.88 | \n5.537493 | \n
| Tlaxcala | \nMan | \n23968.33 | \n1621.2198 | \n21301.47 | \n26635.20 | \n6.764007 | \n
| Tlaxcala | \nWoman | \n22295.82 | \n1473.2391 | \n19872.38 | \n24719.26 | \n6.607692 | \n
| Veracruz De Ignacio De La Llave | \nMan | \n16365.05 | \n1055.8247 | \n14628.24 | \n18101.86 | \n6.451705 | \n
| Veracruz De Ignacio De La Llave | \nWoman | \n14712.45 | \n947.9121 | \n13153.15 | \n16271.74 | \n6.442927 | \n
| Yucatan | \nMan | \n20677.82 | \n1239.0909 | \n18639.55 | \n22716.10 | \n5.992366 | \n
| Yucatan | \nWoman | \n17805.62 | \n1141.8291 | \n15927.34 | \n19683.90 | \n6.412745 | \n
| Zacatecas | \nMan | \n17328.49 | \n1346.8466 | \n15112.96 | \n19544.02 | \n7.772440 | \n
| Zacatecas | \nWoman | \n15690.25 | \n1335.9463 | \n13492.65 | \n17887.85 | \n8.514502 | \n
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Using a Tidy Approach in {R} to Analyze {Mexico's} 2024\n {ENVIPE}},\n date = {2025-01-27},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/estimating-envipe-prevalence/},\n langid = {en}\n}\nThe Standardized Precipitation Index (SPI) quantifies how precipitation deviates from the historical average over multiple time scales (e.g., from 1 to 12 months), allowing for consistent monitoring of droughts and floods across regions and time, making it a valuable tool for early warning and risk assessment systems.
\nIn this tutorial, I will walk through the process of calculating the SPI for Somalia at the subnational level using the {SCI} package in R.
\nThe rainfall dataset used in this tutorial comes from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA) and provides monthly precipitation data for Somalia's subnational districts. It is derived from CHIRPS v2 (Climate Hazards Group InfraRed Precipitation with in situ Station data), using satellite imagery combined with ground-based observations.
\nFor each 10-day period (dekad), the dataset includes the following indicators:
\nrfh: Total rainfall (mm) over 10 daysr1h: 1-month rolling rainfall sum (mm)r3h: 3-month rolling rainfall sum (mm)rfh_avg: Long-term average rainfall for the dekad (mm)r1h_avg: Long-term average for 1-month rolling rainfall (mm)r3h_avg: Long-term average for 3-month rolling rainfall (mm)rfq: Rainfall anomaly (%)r1q: Rainfall anomaly (%)r3q: 3-month rainfall anomaly (%)The data is aggregated at the subnational level and follows administrative boundaries defined by the World Food Programme (WFP), with each unit assigned a unique code (Pcode).
Before we begin, we need to install and load the necessary R packages. We will be using the {tidyverse} package for data wrangling, as well as the {SCI} package for calculating the SPI.
\nlibrary(tidyverse) \nlibrary(SCI) # For Standardized Precipitation Index calculations (climate analysis)\nlibrary(conflicted) # For resolving conflicts between function names in different packages\nlibrary(kableExtra) \n\nlibrary(ggplot2) \nlibrary(ggtext) \nlibrary(showtext) \n\nconflicts_prefer(dplyr::select) # Prefer 'select' function from dplyr package\nconflicts_prefer(dplyr::lag) # Prefer 'lag' function from dplyr package\nconflicts_prefer(dplyr::filter) # Prefer 'filter' function from dplyr packageImportant to use {conflicted} package to specify preferred functions to avoid naming conflicts between {tidyverse} and {SCI} packages.
\nData can be found here in file som-rainfall-adm2-full.csv. For Somalia, data is available from January 1981.
# Read the rainfall data for Somalia`\n# You can also download the data from the HDX platform and read it locally\nrainfall_data <- read.csv(\"https://data.humdata.org/dataset/ed6e1b4b-8094-47e6-bdf7-f6d56fa7abb9/resource/8b333d58-d69e-418c-b5e3-dd86f12eee05/download/som-rainfall-subnat-full.csv\")We preprocess the data to extract relevant information to calculate the SPI. We extract the year and month from the date column, calculate the rainfall anomaly, and group the data by year, month, and district_code.
\nrainfall_data <- rainfall_data |>\n slice(-1) |> # Remove the first row (metadata)\n filter(adm_level == 2) |>\n mutate(\n year = year(date), # Extract year from the 'date' column\n month = month(date), # Extract month from the 'date' column\n district_code = PCODE # Rename the 'ADM2_PCODE' column to 'district_code'\n ) |>\n filter(year < 2025) |> # Filter data up to 2024\n rename(\n rainfall_monthly_total = rfh, # Rename the column 'rfh' to 'rainfall_monthly_total'\n rainfall_monthly_lt_average = rfh_avg # Rename the column 'rfh_avg' to 'rainfall_monthly_lt_average'\n ) |>\n mutate(\n rainfall_monthly_total = as.numeric(rainfall_monthly_total), # Convert rainfall total to numeric\n rainfall_monthly_lt_average = as.numeric(rainfall_monthly_lt_average) # Convert long-term average to numeric\n ) |>\n select(year, month, district_code, rainfall_monthly_total, rainfall_monthly_lt_average) # Select relevant columns\n\n# Group by year, month, and district_code, and calculate mean rainfall values by district\nrainfall_data <- rainfall_data |>\n group_by(year, month, district_code) |>\n summarise(\n rainfall_monthly_total = mean(rainfall_monthly_total, na.rm = TRUE), # Calculate the mean of rainfall total\n rainfall_monthly_lt_average = mean(rainfall_monthly_lt_average, na.rm = TRUE) # Calculate the mean of long-term average\n ) |>\n ungroup() |>\n # Optionally, you can calculate rainfall anomaly as a percentage of the long-term average\n # rainfall_anomaly = 1, indicates no anomaly\n # < 1 indicates below average\n # > 1 indicates above average\n mutate(\n rainfall_anomaly = rainfall_monthly_total / rainfall_monthly_lt_average\n )rainfall_data <- rainfall_data |>\n group_by(district_code) |> # Group by district_code\n mutate(\n spi1 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1, # Start month for SPI calculation\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1, # Start month for SPI calculation\n time.scale = 1, # Time scale for SPI calculation\n distr = \"gamma\", # Distribution for fitting\n p0 = TRUE # Use initial parameter estimates\n )\n ),\n spi3 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1,\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1,\n time.scale = 3,\n distr = \"gamma\",\n p0 = TRUE\n )\n ),\n spi6 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1,\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1,\n time.scale = 6,\n distr = \"gamma\",\n p0 = TRUE\n )\n ),\n spi12 = transformSCI(\n rainfall_monthly_total,\n first.mon = 1,\n obj = fitSCI(rainfall_monthly_total,\n first.mon = 1,\n time.scale = 12,\n distr = \"gamma\",\n p0 = TRUE\n )\n )\n ) |>\n ungroup() |>\n arrange(district_code)The output table looks like this:
\n| year | \nmonth | \ndistrict_code | \nrainfall_monthly_total | \nrainfall_monthly_lt_average | \nrainfall_anomaly | \nspi1 | \nspi3 | \nspi6 | \nspi12 | \n
|---|---|---|---|---|---|---|---|---|---|
| 1981 | \n1 | \nSO1101 | \n0.9417476 | \n1.744876 | \n0.5397218 | \n-1.1043463 | \nNA | \nNA | \nNA | \n
| 1981 | \n2 | \nSO1101 | \n2.0194174 | \n2.027616 | \n0.9959566 | \n0.1537152 | \nNA | \nNA | \nNA | \n
| 1981 | \n3 | \nSO1101 | \n32.8349507 | \n9.696980 | \n3.3861008 | \n2.0730688 | \n2.0600381 | \nNA | \nNA | \n
| 1981 | \n4 | \nSO1101 | \n18.3139157 | \n24.722222 | \n0.7407876 | \n-0.3498783 | \n0.9740066 | \nNA | \nNA | \n
| 1981 | \n5 | \nSO1101 | \n6.7928803 | \n17.754692 | \n0.3825963 | \n-1.1682128 | \n0.3482195 | \nNA | \nNA | \n
| 1981 | \n6 | \nSO1101 | \n2.7443365 | \n4.763215 | \n0.5761522 | \n-1.1812731 | \n-1.5032149 | \n0.1598925 | \nNA | \n
| 1981 | \n7 | \nSO1101 | \n17.3592230 | \n14.047789 | \n1.2357264 | \n0.7686934 | \n-0.9048239 | \n0.3881687 | \nNA | \n
| 1981 | \n8 | \nSO1101 | \n34.8381870 | \n26.758146 | \n1.3019657 | \n0.8281462 | \n0.7648482 | \n0.7461252 | \nNA | \n
| 1981 | \n9 | \nSO1101 | \n29.1682857 | \n19.003020 | \n1.5349290 | \n1.2521098 | \n1.2306838 | \n0.0193027 | \nNA | \n
| 1981 | \n10 | \nSO1101 | \n4.1974111 | \n8.226213 | \n0.5102483 | \n-0.3816705 | \n0.7653264 | \n0.1559530 | \nNA | \n
Drought events are indicated by negative SPI values, with more negative values indicating more severe drought conditions. Conversely, positive SPI values indicate wetter-than-average conditions. For demonstrating purposes, let's visualize the SPI values for Laasqoray district in Sanaag region. One efficient way to visualize SPI values is by using a box plot with SPI categories.
\nfont_add_google(\"Roboto Condensed\", \"Roboto Condensed\")\n\nshowtext_auto()\n\n# Custom theme for the chart\ntheme_boxplot_chart <- function() {\n theme_minimal(\n base_family = \"Roboto Condensed\" \n ) +\n # Custom theme settings\n theme(\n # Axis settings\n axis.title = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 14\n ),\n axis.text = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 12\n ),\n\n # Title settings\n plot.title.position = \"plot\", \n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 24,\n margin = margin(5, 0, 5, 0), \n width = unit(1, \"npc\") \n ),\n\n # Legend settings\n legend.position = \"top\",\n legend.title = element_blank(),\n legend.key.height = unit(0.5, \"cm\"),\n legend.key.width = unit(0.5, \"cm\"),\n legend.spacing.x = unit(0.2, \"cm\"),\n legend.key.spacing = unit(0.5, \"cm\"), \n legend.text = element_text(\n margin = margin(5, 0, 5, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 10\n ),\n legend.direction = \"horizontal\",\n legend.byrow = FALSE,\n\n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n face = \"italic\",\n size = 12,\n hjust = 0\n ),\n plot.background = element_rect(\n color = \"white\",\n fill = \"white\"\n ),\n plot.margin = margin(20, 40, 20, 40)\n )\n}\n\ntitle_chart <- \"1-month Standardized Precipitation Index (SPI-1) in Laasqoray, Somalia for the period 1981-2024\"\ncaption_chart <- paste0(\"**Data:** Own estimations using Somalia rainfall Indicators at Subnational Level (HDX)\",\n \"<br>\", \n \"**Graphic:** Juan Torres Mungu\u00eda\")\n\n#\nspi_colors <- data.frame(\n ymin = c(-5, -2, -1.5, -1, 0, 1, 1.5, 2),\n ymax = c(-2, -1.5, -1, 0, 1, 1.5, 2, 5),\n label = factor(\n c(\n \"Extremely dry\", \"Severely dry\", \"Moderately dry\", \"Near normal\",\n \"Near normal\", \"Moderately wet\", \"Severely wet\", \"Extremely wet\"\n ),\n levels = c(\n \"Extremely dry\", \"Severely dry\", \"Moderately dry\", \"Near normal\",\n \"Moderately wet\", \"Severely wet\", \"Extremely wet\"\n )\n ),\n fill = c(\n \"#F76D5E\", \"#FFAD72\", \"#FFE099\", \"#99FF99\",\n \"#99FF99\", \"#99EAFF\", \"#75D3FF\", \"#3D87FF\"\n )\n)\n\nrainfall_data |>\n filter(district_code == \"SO1503\") |>\n mutate(month = factor(month, levels = 1:12, labels = month.name)) |> # Convert month to factor\n ggplot(aes(x = month, y = spi1)) +\n\n # Add box plot with SPI categories\n geom_rect(\n data = spi_colors, aes(ymin = ymin, ymax = ymax, xmin = -Inf, xmax = Inf, fill = label),\n alpha = 0.2, inherit.aes = FALSE\n ) +\n geom_boxplot(size = 1, outlier.size = 2, outlier.shape = 16) +\n scale_fill_manual(values = setNames(spi_colors$fill, spi_colors$label)) +\n labs(\n title = title_chart,\n caption = caption_chart,\n x = \"\",\n y = \"SPI-1\",\n fill = \"\"\n ) +\n guides(fill = guide_legend(nrow = 1)) + # Number of rows in the legend\n theme_boxplot_chart()showtext_opts(dpi = 320) # Set the resolution of the image 320 dpi is for high-quality images (\"retina\")\nggsave(\n \"spi-somalia.png\",\n dpi = 320,\n width = 12,\n height = 9,\n units = \"in\"\n)\nshowtext_auto(FALSE) # Turn off the showtext functionality
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {How to Compute and Visualize the {Standardized}\n {Precipitation} {Index} {(SPI)} for Drought Analysis},\n date = {2025-02-08},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/spi_somalia/},\n langid = {en}\n}\nPreviously, in Part 1 of this tutorial, we calculated the Standardized Precipitation Index (SPI) for different time scales using the {SCI} package in R. Now, we will map the SPI values for the districts of Somalia using {sf} and {ggplot2} packages in R.
\nThe information of the Somalian sub national administrative boundaries is obtained from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA). The information comes as a shapefile .shp and is available here in a zip file named som_adm_ocha_20250108_AB_SHP.zip. The shapefile contains the administrative boundaries of Somalia at the subnational level. In particular, we will use file som_admbnda_ADM2_ocha_20250108.shp, that contains the shapefile of Somalia at the district level. I downloaded this file and located it in /shp folder.
Before we begin, we need to install and load the necessary R packages.
\n\nData can be found here in file som-rainfall-adm2-full.csv. For Somalia, data is available from January 1981.
I will add the SPI values previously calculated to the shapefile data frame. I saved the SPI values in a RData file named rainfall_data.RData.
load(\"rainfall_data.RData\")\n\n# First, we need to filter the SPI values for the last available month, i.e. December 2024\n# to have only one observation per district\nrainfall_data_dec2024 <- rainfall_data |>\n filter(month == 2 & year == 2024)\n\n# Merge the shapefile with the SPI values\nshp_somalia <- shp_somalia |>\n left_join(rainfall_data_dec2024, by = c(\"district_code\"))font_add_google(\"Roboto Condensed\", \"Roboto Condensed\")\n\nshowtext_auto()\n\ntheme_spi_map <- function() {\n theme_minimal(\n base_family = \"Roboto Condensed\" \n ) +\n # Custom theme settings\n theme(\n\n # Title settings\n plot.title.position = \"plot\", \n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 22,\n margin = margin(5, 0, 5, 0), \n width = unit(1, \"npc\") \n ),\n plot.subtitle = element_textbox(\n color = \"grey50\",\n face = \"bold\",\n size = 14,\n margin = margin(20, 0, 10, 0),\n width = unit(1, \"npc\")\n ),\n\n legend.position = \"top\",\n legend.title = element_blank(),\n legend.key.height = unit(0.1, \"cm\"), \n legend.key.width = unit(0.1, \"cm\"), \n legend.spacing.x = unit(0.1, \"cm\"),\n legend.key.spacing = unit(0.1, \"cm\"), \n legend.text = element_text(\n margin = margin(5, 0, 5, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 10\n ),\n legend.direction = \"horizontal\",\n legend.byrow = FALSE,\n\n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n face = \"italic\",\n size = 10,\n hjust = 0\n ),\n plot.background = element_rect(\n color = \"white\",\n fill = \"white\"\n ),\n plot.margin = margin(20, 40, 20, 40)\n )\n}\n\n# Title, subtitle, and caption for the waffle chart\ntitle_chart <- \"Standardized Precipitation Index (SPI) in Somalia\"\nsubtitle_chart <- \"February 2024\"\ncaption_chart <- \"**Data:** Somalia: Own estimations using Somalia Rainfall Indicators at Subnational Level (HDX) <br> **Graphic:** Juan Torres Mungu\u00eda\"Now, lets create a histogram of the SPI values for the districts of Somalia to use as a legend for the map.
\nspi_colors <- data.frame(\n ymin = c(-5, -2, -1.5, -1, 0, 1, 1.5, 2),\n ymax = c(-2, -1.5, -1, 0, 1, 1.5, 2, 5),\n label = factor(\n c(\n \"Extreme dry\", \"Severe dry\", \"Moderate dry\", \"Near normal\",\n \"Near normal\", \"Moderate wet\", \"Severe wet\", \"Extreme wet\"\n ),\n levels = c(\n \"Extreme dry\", \"Severe dry\", \"Moderate dry\", \"Near normal\",\n \"Moderate wet\", \"Severe wet\", \"Extreme wet\"\n )\n ),\n fill = c(\n \"#F76D5E\", \"#FFAD72\", \"#FFE099\", \"#E5FFE5\",\n \"#E5FFE5\", \"#99EAFF\", \"#75D3FF\", \"#3D87FF\"\n )\n)\n\nshp_somalia <- shp_somalia %>%\n mutate(spi_category = cut(spi1,\n breaks = c(-Inf, spi_colors$ymax),\n labels = spi_colors$label,\n right = TRUE\n )) |> # right = TRUE means that the intervals are closed on the right\n filter(!is.na(spi1))\n\nhistogram_spi <- shp_somalia |>\n filter(!is.na(spi_category)) |>\n ggplot(aes(x = spi_category, fill = spi_category)) +\n geom_bar(color = \"white\") +\n scale_fill_manual(values = setNames(spi_colors$fill, spi_colors$label)) +\n labs(\n title = \"\",\n x = \"\",\n y = \"Number of districts\"\n ) +\n theme_spi_map() +\n theme(\n legend.position = \"none\",\n axis.title = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 7\n ),\n axis.text = element_text(\n color = \"grey10\",\n face = \"bold\",\n size = 7\n ),\n axis.line = element_blank(),\n panel.grid = element_blank()\n ) +\n coord_flip()# Read the shapefile of Somalia with the interior region borders\nshp_somalia_region <- read_sf(\"shp/som_admbnda_adm1_ocha_20250108.shp\")\n\nshp_somalia_region <- st_transform(shp_somalia_region, 3857) # Web Mercator projection\nshp_somalia <- st_transform(shp_somalia, 3857)\n\nmap_spi <- ggplot() +\n geom_sf(\n data = shp_somalia, aes(fill = spi_category), color = \"grey\", linewidth = 0.5\n ) +\n scale_fill_manual(values = setNames(spi_colors$fill, spi_colors$label)) +\n geom_sf(data = shp_somalia_region, fill = NA, color = \"black\", linewidth = 0.9) +\n labs(\n title = title_chart,\n subtitle = subtitle_chart,\n caption = caption_chart,\n x = \"\",\n y = \"\",\n fill = \"\"\n ) +\n theme_spi_map() +\n theme(\n legend.position = \"none\",\n # Axis settings\n axis.title = element_blank(),\n axis.line = element_blank(),\n axis.text = element_blank(),\n axis.ticks = element_blank(),\n panel.grid = element_blank()\n )map_spi +\n inset_element(histogram_spi, left = 0.4, bottom = 0, right = 1, top = 0.4, on_top = FALSE)showtext_opts(dpi = 320) \nggsave(\n \"spi-somalia-map.png\",\n dpi = 320,\n width = 12,\n height = 9,\n units = \"in\"\n)\nshowtext_auto(FALSE)
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {How to Compute and Visualize the {Standardized}\n {Precipitation} {Index} {(SPI)} for Drought Analysis},\n date = {2025-02-14},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/spi_somalia_visualization/},\n langid = {en}\n}\nACLED (Armed Conflict Location & Event Data Project) provides detailed data on political violence and protest events across the globe. This blog post explores the visualization of armed conflict in the Lake Chad Basin, focusing on the dynamics of violence against civilians in the region. Using ACLED's data in combination with shapefiles from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA), we leverage {ggplot2}, {sf} and the {tidyverse} to create insightful maps that highlight the impact of armed conflict.
\nLoad necessary libraries for data processing and visualization.
\nlibrary(tidyverse) # For data manipulation and visualization\nlibrary(sf) # For handling spatial data\nlibrary(ggplot2) # For data visualization\nlibrary(patchwork) # For combining plots\nlibrary(ggtext) # Use rich text elements in ggplot\nlibrary(showtext) # Add custom fonts to plots\nlibrary(ggspatial) # For interacting with spatial data using ggplot2\nlibrary(rnaturalearth) # To provide map data that can be visualized using other R packages\nlibrary(rnaturalearthdata) # To store data for the rnaturalearth package\n\n# Set environment to English (useful for date/time formatting)\ninvisible(Sys.setenv(LANG = \"en\"))\ninvisible(Sys.setlocale(\"LC_TIME\", \"English\"))Accessing ACLED data requires registering on their ACLED Access Portal to generate a unique access key. Once registered, you can download the data using their Data Export Tool. I downloaded all the events occurred in Africa since 1997 and save them in file acled_africa_1997-2025.csv. To have a better understanding of their data, they provide a codebook.
conflict_data <- read.csv(\"acled_africa_1997-2025.csv\")Information of the national and sub national administrative boundaries in the Lake Chad Basin is obtained from the HDX. The information comes as a shapefile .shp and is available here in a zip file named lcb_admbnda_adm2_ocha.zip. The shapefile contains the administrative boundaries of Lake Chad Basin region at the subnational level. In particular, we will use file lcb_admbnda_adm2_ocha.shp. I also downloaded the administrative boundaries of the West and Central Africa region available here.
# Lake Chad Basin Area - Administrative boundaries for Nigeria, Chad, Cameroon and Niger\nshp_lake_chad_basin <- read_sf(\"lcb_admbnda_adm2_ocha/lcb_admbnda_adm2_ocha.shp\")\n\n# West and Central Africa - Administrative boundaries levels 0 - 2 and Settlements\nwest_central_africa <- read_sf(\"wca_admbnda_adm2_ocha/wca_admbnda_adm2_ocha.shp\")As indicated in the ACLED's codebook, geographic location of the events is specified in columns latitude and longitude in the conflict_data object and they are using the EPSG:4326 coordinate reference system. We have to ensure, that both shp_lake_chad_basin and west_central_africa have the same coordinate reference system. To achieve this, we use st_transform() from {sf} package.
# st_transform to convert coordinates to ACLED's data coordinate reference system\nshp_lake_chad_basin <- st_transform(shp_lake_chad_basin,\n crs = 4326\n)\nwest_central_africa <- st_transform(west_central_africa,\n crs = 4326\n)Now, we subset the events to exclusively account for those occurred between 2020 and 2024 in the Lake Chad Basin. To achieve this, we use the {tidyverse} package for data wrangling and {sf} for spatial filtering. By applying the st_filter function, we retain only events located within the Lake Chad Basin area.
conflict_data <- conflict_data |>\n # Keep only events from 2020 to 2024\n filter(year >= 2020 & year < 2025) |>\n # Select only \"Violence against civilians\" events of type \"Attack\"\n filter(event_type == \"Violence against civilians\" & sub_event_type == \"Attack\") |>\n # Keep events from Nigeria, Niger, Cameroon, and Chad\n filter(country %in% c(\"Nigeria\", \"Niger\", \"Cameroon\", \"Chad\")) |>\n # Extract the month name from event_date\n mutate(\n month = month(dmy(event_date),\n label = TRUE\n )\n ) |>\n # Select only relevant columns\n select(\n event_id_cnty, event_date, month, year, event_type, sub_event_type,\n actor1, actor2, location, latitude, longitude, fatalities, notes\n )\n\n# Convert the data to a spatial object using latitude and longitude (WGS 84 CRS)\nconflict_data <- st_as_sf(conflict_data,\n coords = c(\"longitude\", \"latitude\"),\n crs = 4326\n)\n\n# Filter events that fall within the Lake Chad Basin shapefile\nconflict_data_lake_chad <- conflict_data |>\n st_filter(shp_lake_chad_basin)We then processes administrative boundary data to extract relevant regions in West and Central Africa, focusing on the Lake Chad Basin area.
\n# Filter to include only countries in the Lake Chad Basin region\n# I included Central African Republic for visualization purposes\nshp_countries <- west_central_africa |>\n filter(admin0Name %in% c(\"Nigeria\", \"Niger\", \"Chad\", \"Cameroon\", \"Central African Republic\"))\n\n# Create country codes and format administrative division codes\nshp_countries <- shp_countries |>\n mutate(\n country_code = case_when(\n admin0Name == \"Nigeria\" ~ \"NGA\",\n admin0Name == \"Niger\" ~ \"NER\",\n admin0Name == \"Chad\" ~ \"TCD\",\n admin0Name == \"Cameroon\" ~ \"CMR\"\n ),\n # Extract and format adm1 and adm2 codes\n adm1_code = str_pad(str_extract(admin1Pcod, \"\\\\d+\"), 3, pad = \"0\"),\n adm2_code = str_pad(str_sub(admin2Pcod, 5, 6), 3, pad = \"0\"),\n # Create a unique administrative code by combining country and division codes\n adm_code = paste0(country_code, adm1_code, adm2_code)\n ) |>\n # Rename columns for clarity\n rename(\n country = admin0Name,\n adm1 = admin1Name,\n adm2 = admin2Name\n ) |>\n # Select only relevant columns\n select(adm_code, country, adm1, adm2, geometry)\n\n# Identify administrative units within the Lake Chad Basin\nshp_countries <- shp_countries |>\n mutate(lake_chad_basin = ifelse(adm_code %in% shp_lake_chad_basin$Rowcacode2,\n \"Yes\", \"No\"\n ))\n\n# Subset to retain only the administrative units inside the Lake Chad Basin\nshp_lake_chad <- shp_countries |>\n filter(lake_chad_basin == \"Yes\")\n\n# Compute the bounding box for visualization\nbbox <- st_bbox(shp_lake_chad)\nx_min <- bbox[\"xmin\"]\nx_max <- bbox[\"xmax\"]\ny_min <- bbox[\"ymin\"]\ny_max <- bbox[\"ymax\"]\n\n# Aggregate administrative level 1 (regions/states) within the Lake Chad Basin\nshp_adm1_aggregated <- shp_countries |>\n group_by(adm1, lake_chad_basin) |>\n summarise(geometry = st_union(geometry))\n\n# Keep only regions in the Lake Chad Basin and clean names\nshp_adm1_aggregated <- shp_adm1_aggregated |>\n filter(lake_chad_basin == \"Yes\") |>\n mutate(adm1 = str_replace_all(adm1, \"_\", \" \")) |>\n # Compute centroid coordinates for labeling\n mutate(\n x = st_coordinates(st_centroid(geometry))[, 1],\n y = st_coordinates(st_centroid(geometry))[, 2]\n )Now, we create a hexagonal grid over the Lake Chad region with a cell size of approximately 10.69 km and fatalities and attack counts are aggregated per grid cell.
\n# Reproject Lake Chad Basin shapefile to UTM Zone 33N (EPSG:32633)\n# This is important to set a cell size in km\nshp_lake_chad <- st_transform(shp_lake_chad, crs = 32633)\n\n# Create a hexagonal grid over the Lake Chad Basin with a cell size of ~10.69 km\nnc_grid <- st_make_grid(shp_lake_chad,\n cellsize = 10690,\n square = FALSE\n)\n\n# Filter the grid to keep only hexagons intersecting the Lake Chad Basin\nnc_grid_filter <- nc_grid |>\n st_sf() |>\n st_filter(shp_lake_chad)\n\n# Reproject the grid and Lake Chad Basin back to WGS84 (EPSG:4326)\nnc_grid_filter <- st_transform(nc_grid_filter, crs = 4326)\nshp_lake_chad <- st_transform(shp_lake_chad, crs = 4326)\n\n# Spatially join conflict data with the hexagonal grid\nconflict_data_grid <- st_join(nc_grid_filter, conflict_data_lake_chad, left = FALSE)\n\n# Aggregate conflict data by grid cell, summing fatalities and counting attacks\nfatalities_per_grid <- conflict_data_grid |>\n group_by(nc_grid) |>\n summarise(\n fatalities = sum(fatalities, na.rm = TRUE),\n num_attacks = n()\n ) |>\n st_as_sf() |>\n # Compute centroid coordinates for visualization\n mutate(\n x = st_coordinates(st_centroid(nc_grid))[, 1],\n y = st_coordinates(st_centroid(nc_grid))[, 2]\n )We need to relocate the coordinates of the labels for specific regions for better readability.
\n# Adjust label positions for specific administrative regions to improve readability\nshp_adm1_aggregated <- shp_adm1_aggregated |>\n mutate(\n x = case_when(\n adm1 == \"Extreme Nord\" ~ x - 0.1, # Shift label slightly left\n TRUE ~ x\n ),\n y = case_when(\n adm1 == \"Extreme Nord\" ~ y - 0.5, # Move label downward\n adm1 == \"Adamawa\" ~ y + 0.3, # Move label upward\n adm1 == \"Borno\" ~ y - 0.4, # Adjust label downward\n TRUE ~ y\n )\n )\n\n# Aggregate country geometries by merging administrative boundaries\n# Dissolve some shapes\nshp_countries_countries <- shp_countries |>\n group_by(country) |>\n summarise(geometry = st_union(geometry)) |>\n st_cast(\"POLYGON\") |>\n group_by(country) |>\n # Retain the largest polygon for each country\n filter(st_area(geometry) == max(st_area(geometry))) \n\n# Define manual label positions for country names\ncountry_labels <- tibble(\n country = c(\"Nigeria\", \"Niger\", \"Chad\", \"Cameroon\"),\n x = c(10.15, 11, 15, 13.6), # X-coordinates for label placement\n y = c(10, 15, 14.2, 9) # Y-coordinates for label placement\n)We first create an empty map of the countries and highlight the Lake Chad Basin area using {ggplot2}.
\n# Automatically use showtext for handling fonts\nshowtext_auto()\n\n# Add Google font \"Voltaire\"\nfont_add_google(\"Voltaire\", \"Voltaire\")\n\nlake_chad_map <- ggplot() +\n\n # Plot country borders with a light fill color\n geom_sf(\n data = shp_countries,\n fill = \"#ECEFF1\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n ) +\n\n # Plot Lake Chad Basin area with a darker fill color\n geom_sf(\n data = shp_lake_chad,\n fill = \"#E0E0E0\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n )Now, we add the fatalities information.
\nlake_chad_map <- lake_chad_map +\n # Plot fatalities per grid with a color gradient\n geom_sf(\n data = fatalities_per_grid,\n aes(fill = fatalities),\n color = NA,\n alpha = 1\n ) +\n scale_fill_gradient(name = \"Fatalities\", \n low = \"#FCC5C0\", high = \"#67001F\")We add different color to subnational and national boundaries and add their corresponding labels.
\nlake_chad_map <- lake_chad_map +\n # Plot aggregated administrative region boundaries with white lines\n geom_sf(\n data = shp_adm1_aggregated,\n fill = NA,\n color = \"white\",\n linewidth = 0.5\n ) +\n\n # Plot country boundaries with black lines\n geom_sf(\n data = shp_countries_countries,\n fill = NA,\n color = \"black\",\n linewidth = 1.25\n ) +\n # Add labels for administrative regions (e.g., \"Extreme Nord\")\n geom_text(\n data = shp_adm1_aggregated,\n aes(\n x = x, y = y,\n # Wrap text if needed for long region names\n label = str_wrap(adm1, 1) \n ),\n colour = \"#999\",\n fontface = \"bold\",\n size = 3,\n family = \"Voltaire\"\n ) +\n \n # Add labels for countries with a larger font\n geom_label(\n data = country_labels,\n aes(x = x, y = y, label = country),\n colour = \"black\",\n fontface = \"bold\",\n size = 4,\n family = \"Voltaire\"\n )Finally, we center the map at the Lake Chad Basin area and add format for a better visualization.
\n# Title, subtitle, and caption for the waffle chart\ntitle_chart <- \"Armed conflict in the Lake Chad Basin, 2020-2024\"\nsubtitle_chart <- \"Fatalities in attacks perpetrated by organized armed actors against civilians\"\n\nlake_chad_map <- lake_chad_map +\n\n labs(\n title = title_chart,\n subtitle = subtitle_chart,\n caption = \"**Data:** Armed Conflict Location & Event Data (ACLED) <br> **Graphic:** Juan Torres Mungu\u00eda\",\n x = \"\",\n y = \"\",\n fill = \"\"\n ) +\n\n # Set the map's coordinate limits and prevent expansion\n coord_sf(\n xlim = c(x_min, x_max),\n ylim = c(y_min, y_max),\n expand = FALSE\n ) +\n\n # Add a scale bar in the bottom left\n annotation_scale(location = \"bl\", width_hint = 0.2) +\n\n # Add a north arrow in the top left corner\n annotation_north_arrow(\n location = \"tl\",\n which_north = \"true\",\n style = north_arrow_fancy_orienteering,\n height = unit(1, \"cm\"),\n width = unit(1, \"cm\")\n ) +\n\n # Customize the legend\n guides(\n fill = guide_legend(position = \"bottom\")\n ) +\n\n # Set the base theme without any background elements\n theme_void(base_family = \"Voltaire\") +\n\n # Customize the map's text and legend styling\n theme(\n # Title settings\n plot.title.position = \"plot\", # Position of the title\n plot.title = element_textbox(\n color = \"black\",\n face = \"bold\",\n size = 22,\n margin = margin(5, 3, 5, 3), # Top, right, bottom, left\n # Width of the title, npc == 1 corresponds to the full width of the plot\n width = unit(1, \"npc\") \n ),\n plot.subtitle = element_textbox(\n color = \"grey50\",\n face = \"bold\",\n size = 14,\n margin = margin(20, 3, 10, 3),\n width = unit(1, \"npc\")\n ),\n\n # Caption settings\n plot.caption = element_markdown(\n color = \"grey70\",\n size = 10,\n margin = margin(10, 3, 5, 3)\n ),\n legend.title = element_text(hjust = 0.5, color = \"black\", face = \"bold\"),\n legend.text = element_text(color = \"black\", face = \"bold\"),\n legend.direction = \"horizontal\",\n legend.title.position = \"top\",\n legend.text.position = \"bottom\",\n legend.key.size = unit(0.65, units = \"cm\"),\n legend.background = element_rect(fill = NA, color = NA),\n legend.margin = margin(t = 5, r = 10, b = 5, l = 10),\n plot.background = element_rect(fill = \"#ECEFF1\", color = NA)\n )We create a zoomed-in map of the Lake Chad Basin inset within the African continent, using ggplot2. It includes the ne_countries function from the {rnaturalearth} package to retrieve the map of Africa. The map of the Lake Chad Basin is overlaid on top of Africa, and a rectangle is drawn to highlight the region of interest (Lake Chad Basin). The map is customized with a white background and black borders.
\n# Fetch African countries map from the rnaturalearth package\nafrica <- ne_countries(continent = \"Africa\", returnclass = \"sf\")\n\n# Create a zoomed-in map of Lake Chad Basin within the African map\nlake_chad_map_zoom <- ggplot() +\n\n # Plot the Africa map with a light fill color\n geom_sf(\n data = africa,\n fill = \"#ECEFF1\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n ) +\n\n # Plot Lake Chad Basin area with a darker fill color\n geom_sf(\n data = shp_lake_chad,\n fill = \"#E0E0E0\",\n color = NA,\n linewidth = 0,\n inherit.aes = FALSE\n ) +\n\n # Add a rectangle highlighting the Lake Chad Basin on the Africa map\n geom_rect(\n aes(xmin = x_min, xmax = x_max, ymin = y_min, ymax = y_max),\n color = \"black\",\n fill = NA,\n linewidth = 0.2\n ) +\n\n # Remove grid and axes for a clean map appearance\n theme_void() +\n\n # Customize borders and background of the map\n theme(\n panel.border = element_rect(\n color = \"black\",\n fill = NA,\n size = 0.8\n ),\n plot.background = element_rect(fill = \"white\")\n )Finally, we combine the main map of the Lake Chad Basin with an inset map of Africa showing the location of the Lake Chad Basin. The inset map is positioned within the main map using the inset_element() function of the {patchwork} package. Additionally, we set the resolution of the saved image to 320 dpi (high resolution) for better quality. Finally, the map is saved as a PNG file using ggsave()
# locate main and inset map\nlake_chad_map +\n inset_element(\n align_to = \"plot\",\n lake_chad_map_zoom,\n right = 1,\n left = 0.84,\n top = 1,\n bottom = 0.91,\n on_top = TRUE\n )# Set the resolution of the image 320 dpi is for high-quality images (\"retina\")\nshowtext_opts(dpi = 320) \nggsave(\n \"lake-chad-map.png\",\n dpi = 320,\n width = 5,\n height = 12,\n units = \"in\"\n)\nshowtext_auto(FALSE)
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Visualizing Armed Conflict Using {ACLED's} Data and {HDX}\n Shapefiles},\n date = {2025-03-09},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/visualizing-conflict-lake-chad-basin/},\n langid = {en}\n}\nWe begin by loading the necessary R packages:
\n\nWe retrieve the latest dataset of disease outbreaks from GitHub:
\n# Define the API URL from GitHub containing the latest dataset update\nurl_api <- \"https://api.github.com/repos/jatorresmunguia/disease_outbreak_news/contents/Last%20update\"\n\n# Retrieve the list of files in the \"Last update\" folder\nfile_list <- fromJSON(content(GET(url_api), as = \"text\"))$name\n\n# Filter the most recent file that starts with \"outbreaks\"\nlast_file <- file_list[grepl(file_list, pattern = paste0(\"^outbreaks\"))]\n\n# Select the file with a .csv extension\nrdata_file <- last_file[grepl(\".csv$\", last_file)]\n\n# Construct the CSV file URL from GitHub and read the data into a dataframe\nlast_version <- read.csv(paste0(\"https://raw.githubusercontent.com/jatorresmunguia/disease_outbreak_news/refs/heads/main/Last%20update/\", rdata_file), header = TRUE)\n\n# Remove unnecessary ICD-11 classification columns\nlast_version <- last_version |> select(-c(icd11c1, icd11c2, icd11c3, icd11l1, icd11l2, icd11l3))We integrate additional geographic classifications from WHO and UNSD.
\n# Define the URL for the WHO regional classification dataset\nurl_who_class <- \"https://cdn.who.int/media/docs/default-source/air-pollution-documents/air-quality-and-health/un-agencies-region-classification-for-country.xlsx?sfvrsn=289af35f_3\"\n\n# Create a temporary file to store the downloaded Excel file\ntemp_file <- tempfile(fileext = \".xlsx\")\n\n# Download the WHO classification file and save it to the temporary file\nGET(url_who_class, write_disk(temp_file, overwrite = TRUE))Response [https://cdn.who.int/media/docs/default-source/air-pollution-documents/air-quality-and-health/un-agencies-region-classification-for-country.xlsx?sfvrsn=289af35f_3]\n Date: 2025-04-28 06:21\n Status: 200\n Content-Type: application/vnd.openxmlformats-officedocument.spreadsheetml.sheet\n Size: 28.9 kB\n<ON DISK> C:\\Users\\jator\\AppData\\Local\\Temp\\RtmpiGD1DX\\file990423923f25.xlsx# Read the WHO classification data from the downloaded Excel file\nwho_class <- read_excel(temp_file)\n\n# Select relevant columns (ISO country code and WHO region name)\n# Rename columns for consistency with the main dataset\nlast_version <- who_class |>\n select(`ISO Country code`, `WHO Region name2`) |>\n rename(iso3 = `ISO Country code`, who_region = `WHO Region name2`) |>\n right_join(last_version, by = \"iso3\") # Merge with the main dataset based on the ISO3 country code# Read the UNSD regional classification dataset from an Excel file\n# The dataset is available at: https://unstats.un.org/unsd/methodology/m49/\nunsd_class <- read_excel(path = \"UNSD \u2014 Methodology.xlsx\")\n\n# Select relevant columns: ISO-alpha3 country code, region name, and sub-region name\n# Rename columns to match the main dataset\nlast_version <- unsd_class |>\n select(`ISO-alpha3 Code`, `Region Name`, `Sub-region Name`) |>\n rename(\n iso3 = `ISO-alpha3 Code`,\n unsd_region = `Region Name`,\n unsd_subregion = `Sub-region Name`\n ) |>\n right_join(last_version, by = \"iso3\") # Merge with the main dataset based on the ISO3 country codeTo facilitate better organization, we:
\n# Arrange the dataset in descending order by Year, then by ISO3 country code, and ICD-10 4-character code\nlast_version <- last_version |>\n arrange(desc(Year), iso3, icd104c) |>\n # Create a unique outbreak ID by concatenating Year, ISO3 country code, and ICD-10 4-character code\n mutate(id_outbreak = paste0(Year, iso3, icd104c)) |>\n # Add an empty row at the beginning of the dataset (for later HXL hashtag insertion)\n add_row(.before = 1)We define HXL hashtags for each relevant column:
\n# Define the mapping of column names to their corresponding HXL hashtags\ncols2hxl <- c(\n id_outbreak = \"#Year+iso3+icd4\", # Unique outbreak ID\n Country = \"#country+name\", # Country name\n iso2 = \"#country+code+iso2\", # ISO 2-letter country code\n iso3 = \"#country+code+iso3\", # ISO 3-letter country code\n unsd_region = \"#region+unsd\", # UNSD region\n unsd_subregion = \"#subregion+unsd\", # UNSD subregion\n who_region = \"#region+who\", # WHO region\n Year = \"#date+year\", # Year of the outbreak\n icd10n = \"#disease+name+icd\", # Disease name (ICD-10)\n icd103n = \"#disease+name+icd3\", # Disease name (ICD-10 3-character)\n icd104n = \"#disease+name+icd4\", # Disease name (ICD-10 4-character)\n icd10c = \"#disease+code+icd\", # Disease code (ICD-10)\n icd103c = \"#disease+code+icd3\", # Disease code (ICD-10 3-character)\n icd104c = \"#disease+code+icd4\", # Disease code (ICD-10 4-character)\n Disease = \"#disease+name\", # Disease name\n DONs = \"#news+id\", # News ID\n Definition = \"#x_disease+definition\" # Disease definition\n)\n\n# Add HXL hashtags to the relevant columns in the dataset\nlast_version <- last_version |>\n mutate(across(all_of(names(cols2hxl)), ~ replace(., 1, cols2hxl[cur_column()]))) |>\n # Select the relevant columns for the final dataset\n select(id_outbreak, Year, icd10n, icd103n, icd104n, icd10c, icd103c, icd104c, Disease, Definition, Country, iso2, iso3, unsd_region, unsd_subregion, who_region, DONs)Finally, we save the dataset as an Excel file:
\nwrite_xlsx(x = list(Data = last_version), path = \"disease_outbreaks_HDX.xlsx\")@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Working with {Humanitarian} {Exchange} {Language} {(HXL)}\n {Hashtags}},\n date = {2025-03-26},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/hdx-tags-to-outbreaks-data/},\n langid = {en}\n}\nAccording to the 2021-22 Domestic Violence module of the Kingdom of Cambodia's Demographic and Health Survey (CDHS), over 21% of Cambodian women aged 15\u201349 who have ever had a husband (or intimate partner), have experienced physical, sexual, economic, and/or emotional intimate partner violence in their lifetime. This figures varies significantly across the country, with some provinces reporting rates as high as 45% while others report rates as low as 6%.
\nA better understanding of the geographical distribution of intimate partner violence (IPV) in Cambodia can help policymakers and organizations to target their efforts more effectively. In this blog post, I will show how to create choropleth maps using R and shapefiles to visualize the prevalence of IPV in Cambodia at the provincial level. Additionally, I will use the {patchwork} package to combine these maps into the same graphic. The result will look like this:
\n
Information comes from the document \"Further Data Analysis Report. From the Cambodia Demographic and Health Survey 2021-2022\", available here. In particular, I will use Table A.1: Intimate-partner violence by background characteristics.
\nFirst, we install and load the necessary R packages.
\nlibrary(tidyverse) \nlibrary(sf) # For handling shapefiles\nlibrary(readxl) # To read Excel files\nlibrary(kableExtra) # For data tables\nlibrary(ggplot2) \nlibrary(ggtext) \nlibrary(showtext) \nlibrary(patchwork) # To combine multiple ggplots\nlibrary(colorspace) # For color palettes\nlibrary(grid) # To add lines and text elements to graphical objectsI saved the information from Table A.1: Intimate-partner violence by background characteristics in the Excel file intimate-partner-violence-cambodia.xlsx.
\nipv_cambodia <- read_excel(\"intimate-partner-violence-cambodia.xlsx\")This information looks like this:
\nipv_cambodia |>\n kbl(\n caption = \"Lifetime experience of specific acts of intimate partner physical violence committed by their current or most recent partner by provinces, Cambodia 2021-22\"\n ) |>\n kable_paper(\"hover\", full_width = F)| Region | \nEmotional violence | \nPhysical violence | \nSexual violence | \nPhysical or sexual | \nPhysical or sexual or emotional | \nNumber of women who ever had a husband/intimate partner | \n
|---|---|---|---|---|---|---|
| Banteay Meanchey | \n31.3 | \n18.2 | \n4.6 | \n20.2 | \n34.9 | \n235 | \n
| Battambang | \n24.9 | \n16.1 | \n6.3 | \n19.8 | \n30.4 | \n378 | \n
| Kampong Cham | \n7.7 | \n5.2 | \n1.1 | \n5.8 | \n11.5 | \n352 | \n
| Kampong Chhnang | \n19.2 | \n6.7 | \n1.7 | \n6.7 | \n19.2 | \n213 | \n
| Kampong Speu | \n16.6 | \n9.5 | \n2.6 | \n10.1 | \n19.1 | \n362 | \n
| Kampong Thom | \n32.4 | \n12.1 | \n5.0 | \n14.4 | \n35.1 | \n242 | \n
| Kampot | \n20.8 | \n6.4 | \n4.0 | \n9.0 | \n23.0 | \n218 | \n
| Kandal | \n23.9 | \n6.2 | \n2.0 | \n6.5 | \n23.9 | \n423 | \n
| Koh Kong | \n5.4 | \n3.8 | \n0.4 | \n3.8 | \n6.2 | \n41 | \n
| Kratie | \n13.7 | \n9.1 | \n2.8 | \n10.2 | \n15.9 | \n138 | \n
| Mondul Kiri | \n42.7 | \n19.0 | \n3.6 | \n19.4 | \n45.2 | \n34 | \n
| Phnom Penh | \n10.3 | \n5.1 | \n1.3 | \n5.8 | \n11.0 | \n904 | \n
| Preah Vihear | \n38.1 | \n20.2 | \n8.5 | \n22.4 | \n40.0 | \n108 | \n
| Prey Veng | \n25.0 | \n7.2 | \n1.7 | \n8.6 | \n26.3 | \n366 | \n
| Pursat | \n10.8 | \n3.6 | \n0.7 | \n3.6 | \n10.8 | \n111 | \n
| Ratanak Kiri | \n24.7 | \n6.6 | \n0.8 | \n7.1 | \n25.7 | \n96 | \n
| Siemreap | \n21.5 | \n12.0 | \n2.9 | \n13.5 | \n24.6 | \n482 | \n
| Preah Sihanouk | \n10.3 | \n3.4 | \n1.6 | \n4.0 | \n10.8 | \n75 | \n
| Stung Treng | \n29.3 | \n12.5 | \n7.5 | \n16.3 | \n31.5 | \n66 | \n
| Svay Rieng | \n10.9 | \n8.3 | \n2.1 | \n8.8 | \n13.6 | \n217 | \n
| Takeo | \n18.4 | \n7.4 | \n1.2 | \n8.0 | \n19.3 | \n344 | \n
| Otdar Meanchey | \n15.1 | \n6.1 | \n0.0 | \n6.1 | \n15.4 | \n77 | \n
| Kep | \n8.0 | \n7.2 | \n1.2 | \n7.7 | \n11.0 | \n18 | \n
| Pailin | \n11.3 | \n9.5 | \n1.7 | \n9.5 | \n14.2 | \n29 | \n
| Tboung Khmum | \n11.1 | \n7.3 | \n2.3 | \n7.3 | \n14.0 | \n252 | \n
Now, we load geographical data of the administrative division of Cambodia. These data are obtained from the Humanitarian Data Exchange (HDX) initiative of the United Nations Office for the Coordination of Humanitarian Affairs (OCHA). The information comes as a shapefile .shp and is available here in a zip file named som_adm_ocha_20250108_AB_SHP.zip. I downloaded this file and located it in /khmadm_gov_20181004_ab_shp folder.
shp_cambodia <- read_sf(\"khmadm_gov_20181004_ab_shp/khm_admbnda_adm1_gov_20181004.shp\")Next step is to combine the geographic information from the shapefile, shp_cambodia, with the data frame with the information on IPV, ipv_cambodia. To achieve this goal, we need to make sure that both datasets have a common variable to merge them. In this case, we will use the province name ADM1_EN in the shapefile and Region in the data frame. Then, I compare if the names of the provinces in both datasets match:
setdiff(shp_cambodia$ADM1_EN, ipv_cambodia$Region)[1] \"Oddar Meanchey\"The output indicates that there is one discrepancy. In shp_cambodia, the name of the province \"Oddar Meanchey\" corresponds to \"Otdar Meanchey\" in ipv_cambodia. I will rename this province in the shapefile to match the name in the data from the official report.
After this modification, I join both datasets using a left join to keep all the provinces in the shapefile and add the corresponding IPV values from the data frame:
\n\nMy idea is to develop three maps, each representing one type of violence against women and girls: emotional, physical, and sexual. I will first create the map corresponding to emotional violence.
\n# More font alternatives can be found here https://fonts.google.com/\nfont_add_google(\"Scada\", \"Scada\") \n\nshowtext_auto()\n\nmap_ipv_emotional_cambodia <- ggplot() +\n geom_sf(\n data = shp_cambodia, \n aes(fill = `Emotional violence`), \n color = \"white\", # Color of the province borders\n linewidth = 0.5\n ) +\n geom_sf_text(\n data = shp_cambodia,\n aes(label = `Emotional violence`),\n check_overlap = TRUE, # To avoid overlaping text in the plot\n size = 4, \n color = \"black\",\n fontface = \"bold\"\n ) +\n scale_fill_gradient2(\n low = \"#FFFE9EB3\",\n mid = \"#C53270B3\",\n high = \"#040404B3\",\n midpoint = 22.5,\n limits = c(0, 45),\n breaks = c(0, 15, 30, 45),\n labels = c(\"0\", \"15\", \"30\", \"45\")\n ) +\n labs(\n title = \"Emotional\", # Title of this subplot\n subtitle = \"\",\n caption = \"\",\n x = \"\",\n y = \"\",\n fill = \"Percentage of victims (%)\" # Title of the legend\n ) +\n theme_minimal(\n base_family = \"Scada\"\n ) +\n theme(\n plot.margin = margin(0, 40, 0, 40),\n \n plot.background = element_blank(), # No color in the background\n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 22,\n hjust = 0.5,\n margin = margin(5, 0, 5, 0)\n )\n )Now, I replicate the same format for physical and sexual violence.
\nmap_ipv_physical_cambodia <- ggplot() +\n geom_sf(\n data = shp_cambodia, \n aes(fill = `Physical violence`), \n color = \"white\", \n linewidth = 0.5\n ) +\n geom_sf_text(\n data = shp_cambodia,\n aes(label = `Physical violence`),\n check_overlap = TRUE,\n size = 4, \n color = \"black\",\n fontface = \"bold\"\n ) +\n scale_fill_gradient2(\n low = \"#FFFE9EB3\",\n mid = \"#C53270B3\",\n high = \"#040404B3\",\n midpoint = 22.5,\n limits = c(0, 45)\n ) +\n labs(\n title = \"Physical\",\n subtitle = \"\",\n caption = \"\",\n x = \"\",\n y = \"\",\n fill = \"\"\n ) + \n guides(fill = \"none\") +\n theme_minimal(\n base_family = \"Scada\"\n ) +\n theme(\n plot.margin = margin(0, 40, 0, 40),\n \n plot.background = element_blank(),\n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 22,\n hjust = 0.5,\n margin = margin(5, 0, 5, 0)\n )\n )\n\nmap_ipv_sexual_cambodia <- ggplot() +\n geom_sf(\n data = shp_cambodia, \n aes(fill = `Sexual violence`), \n color = \"white\", \n linewidth = 0.5\n ) +\n geom_sf_text(\n data = shp_cambodia,\n aes(label = `Sexual violence`),\n check_overlap = TRUE,\n size = 4, \n color = \"black\",\n fontface = \"bold\"\n ) +\n scale_fill_gradient2(\n low = \"#FFFE9EB3\",\n mid = \"#C53270B3\",\n high = \"#040404B3\",\n midpoint = 22.5,\n limits = c(0, 45)\n ) +\n labs(\n title = \"Sexual\",\n subtitle = \"\",\n caption = \"\",\n x = \"\",\n y = \"\",\n fill = \"\"\n ) + \n guides(fill = \"none\") +\n theme_minimal(\n base_family = \"Scada\"\n ) +\n theme(\n plot.margin = margin(0, 40, 0, 40),\n \n plot.background = element_blank(),\n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 22,\n hjust = 0.5,\n margin = margin(5, 0, 5, 0)\n )\n )Now, using {patchwork} package, I combine the three plots into one single representation. I also use plot_annotation() function to add a title, subtitle, and caption to the overall design.
# Design display two rows, one with the three maps, and the second one only with the common guide\nmap_ipv_cambodia <- (map_ipv_emotional_cambodia + map_ipv_physical_cambodia + map_ipv_sexual_cambodia) / guide_area() + \n # Custom title, subtitle, and caption\n plot_annotation(\n title = \"Intimate partner violence in Cambodia, 2021-22\",\n subtitle = \"Proportion of ever-partnered women and girls aged 15-49 who have been victims of violence perpetrated by their intimate partner, by type of violence.\",\n caption = \"**Data:** 2021-22 Cambodia Demographic and Health Survey (CDHS) <br> **Chart:** Juan Torres Mungu\u00eda\",\n theme = theme(\n # Title settings\n plot.title.position = \"plot\", \n plot.title = element_text(\n color = \"black\",\n face = \"bold\",\n size = 30,\n hjust = 0,\n margin = margin(5, 0, 5, 50)\n ),\n \n plot.subtitle = element_text(\n color = \"grey20\",\n size = 25,\n hjust = 0,\n margin = margin(5, 0, 40, 50)\n ),\n \n # Caption settings\n # I use element_markdown() because the caption has markdown text\n plot.caption = element_markdown(\n color = \"black\",\n size = 15,\n hjust = 0,\n margin = margin(80, 0, 5, 50)\n )\n )) +\n\n plot_layout(guides = \"collect\",\n heights = unit(c(13, 1),\n c('cm', 'cm'))\n ) + \n theme_minimal(\n base_family = \"Scada\"\n ) &\n # Custom theme settings\n theme(\n legend.background = element_blank(),\n legend.position = \"bottom\",\n legend.title = element_text(\n margin = margin(5, 0, 10, 0),\n face = \"bold\",\n color = \"grey10\",\n size = 14,\n hjust = 0.5\n ),\n legend.title.position = \"top\",\n legend.text = element_text(\n margin = margin(5, 5, 5, 5), \n face = \"bold\",\n color = \"grey10\",\n size = 14\n ),\n legend.key.width = unit(1.4, \"cm\"),\n legend.key.height = unit(1, \"cm\"),\n legend.direction = \"horizontal\",\n\n plot.background = element_rect(\n color = \"#E5E5E5\",\n fill = \"#E5E5E5\"),\n # Axis settings\n axis.title = element_blank(),\n axis.line = element_blank(),\n axis.text = element_blank(),\n axis.ticks = element_blank(),\n panel.grid = element_blank()\n \n )Finally, I use grid.lines() from {grid} package to add a horizontal line after the subtitle.
map_ipv_cambodia\n\ngrid.lines(\n x = c(0.04, 0.96), # Values should be between 0 and 1\n y = 0.88, # Value should be between 0 and 1\n gp = gpar(col = \"#C53270B3\", \n lwd = 4)\n)
@online{torres_mungu\u00eda2025,\n author = {Torres Mungu\u00eda, Juan Armando},\n title = {Three Maps to Examine the Prevalence of Intimate Partner\n Violence {(IPV)} in {Cambodia}},\n date = {2025-05-10},\n url = {https://juan-torresmunguia.netlify.app/blog/posts/map-intimate-partner-violence-cambodia/},\n langid = {en}\n}\n