Average color of each country (using satellite imagery)

Couple of weeks ago we have noticed the article (opens new window) getting attention around twitter (opens new window). As it is using our WMS (opens new window) EOxMaps endpoint we have asked Erin Davis (opens new window) if we could host her idea here. At EOX we endorse brilliant and unique ideas that utilize our tools and data. Thank you for this most interesting contribution. Now follows the original article.

---

I wrote that I wanted to expand creatively beyond maps in 2021, but here we are, halfway in, and half my posts are maps.

In my defense, I made these last year and just haven’t gotten around to posting them yet! They’ve been sitting in a folder since December, mocking me. It’s time to get these guys off my computer and out into the world so I can move on to other things.

Since it’s been so long, I really don’t remember what inspired me to make these. I’m sure memories of this amazing NYT Upshot article (opens new window) were involved. The difference here is I didn’t even try to tell a story with these maps… I just made them as something pretty to look at.

Average colors of the USA

A question: do you think it’s best practice to label US states? I usually haven’t in my personal work, but I do in my day job… the labeled look might be growing on me.

Average colors of the world

How I did it

To do this, you’ll need a few things:

1. A shapefile of the areas you want the colors for (from e. g.: https://www.naturalearthdata.com/ (opens new window))
2. QGIS (opens new window) installed (or GDAL/Python, but installing QGIS will set it all up for you)
3. R and RStudio (opens new window) installed

In short, my process was the following:

1. Wrote a script (in R) to find the bounding boxes of all the areas I was interested in
2. In R, downloaded pictures of those areas from Sentinel-2
3. In R, wrote a script that created a series of GDAL (opens new window) commands to:
   1. Georeference the downloaded image to create a geotiff
   2. Crop the geotiff to the borders of the area
   3. Re-project the geotiff to a sensible projection
4. Ran that GDAL script
5. In R, converted the geotiffs back to pngs, and found the average color of the png

Basically, I did everything in R, including creating the GDAL script. I just briefly stepped outside of R to run that script then jumped back in to finish the project. What the GDAL script does is possible within R using the qgisprocess (opens new window) package, but I’m pretty sure (based on 0 testing tbh) that GDAL is faster.

Side note: I feel really embarrassed lifting the curtain and showing my wacky processes to get viz done. I’ve often avoided sharing, especially in professional situations, for fear that better-informed people will laugh at me. But I think it’s good to have more people sharing their code so we can all improve. All my GIS knowledge (and quite a lot of my general coding knowledge) is self-taught, and I owe a huge debt to people online who bothered to explain how they did things. My point? I don’t know… I guess if you know of a better way to do this, please let me know! But also I was really proud I figured out this much on my own.

 

library(magick)
library(sf)
library(dplyr)

1
2
3
4
5
6
7

Get your shapefile loaded and prepped:

setwd("C:/Users/Erin/Documents/DataViz/AvgColor/Temp/")
 

dir.create(file.path(getwd(), "satellite"), showWarnings = FALSE)
dir.create(file.path(getwd(), "shapes"), showWarnings = FALSE)
 


shapes <- read_sf("cb_2018_us_state_20m.shp")
 

for (i in 1:nrow(shapes)) {
  shape <- shapes[i, ]
  write_sf(shape, paste0("./shapes/", shape$NAME, ".shp"), delete_datasource = TRUE)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Download the satellite imagery

shapes <- read_sf("cb_2018_us_state_20m.shp")
 

for (i in 1:nrow(shapes)) {
  shape <- shapes[i, ] %>% st_transform(4326)
  box <- st_bbox(shape)
   
  url <- paste0("https://tiles.maps.eox.at/wms?service=wms&request=getmap&version=1.1.1&layers=s2cloudless-2019&",
                "bbox=", box$xmin-.05, ",", box$ymin-.05, ",", box$xmax+.05, ",", box$ymax+.05,
                "&width=4096&height=4096&srs=epsg:4326")
   
  
  
  download.file(url, paste0("./satellite/", gsub("/", "", shape$NAME), ".jpg"), mode = "wb")
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

Create a GDAL script to transform those square satellite images into cropped, reprojected geotifs.

Note: for some reason, GDAL refuses to deal with any folders other than my temp directory. I couldn’t figure out how to get it to successfully write to my working directory, so I just did stuff in the temp one. Minor annoyance, but what can you do.

if (file.exists("script.txt")) { file.remove("script.txt")}
 

for (i in 1:nrow(shapes)) {
  shape <- shapes[i, ]
  box <- st_bbox(shape)
   
  file_origin <- paste0(getwd(), '/satellite/', shape$NAME, '.jpg')
  file_temp <- paste0(tempdir(), "/", shape$NAME, '.jpg')
  file_final <- paste0(tempdir(), "/", shape$NAME, '_ref.tif')
   
  cmd1 <- paste0("gdal_translate -of GTiff",
                 " -gcp ", 0, " ", 0, " ", box$xmin, " ", box$ymax,
                 " -gcp ", 4096, " ", 0, " ", box$xmax, " ", box$ymax,
                 " -gcp ", 0, " ", 4096, " ", box$xmin, " ", box$ymin,
                 " -gcp ", 4096, " ", 4096, " ", box$xmax, " ", box$ymin,
                 ' "', file_origin, '"',
                 ' "', file_temp, '"')
   
  cmd2 <- paste0("gdalwarp -r near -tps -co COMPRESS=NONE  -t_srs EPSG:4326",
                 ' "', file_temp, '"',
                 ' "', file_final, '"')
   
  write(cmd1,file="script.txt",append=TRUE)
  write(cmd2,file="script.txt",append=TRUE)
}
 

for (i in  1:nrow(shapes)) {
  shape <- shapes[i, ]
  file_shp <- paste0(getwd(), '/shapes/', shape$NAME, '.shp')
  file_orig <-  paste0(tempdir(), "/", shape$NAME, '_ref.tif')
  file_crop <-  paste0(tempdir(), "/", shape$NAME, '_crop.tif')
   
  cmd <- paste0("gdalwarp -cutline",
                ' "', file_shp, '"',
                " -crop_to_cutline -of GTiff  -dstnodata 255",
                ' "', file_orig, '"',
                ' "', file_crop, '"',
                " -co COMPRESS=LZW -co TILED=YES --config GDAL_CACHEMAX 2048 -multi")
   
  write(cmd,file="script.txt",append=TRUE)
}
 

for (i in  1:nrow(shapes)) {
  shape <- shapes[i, ]
   
  center <- st_centroid(shape) %>% st_coordinates
  utm <- as.character(floor((center[1] + 180) / 6) + 1)
  if (nchar(utm) == 1) {utm <- paste0("0", utm)}
  if (center[2] >=0) {
    epsg <- paste0("326", utm)
  } else {
    epsg <- paste0("327", utm)
  }
   
  file_orig <- paste0(tempdir(), "/", shape$NAME, '_crop.tif')
  file_mod <- paste0(tempdir(), "/", shape$NAME, '_reproj.tif')
   
  cmd <- paste0("gdalwarp -t_srs EPSG:", epsg,
                ' "', file_orig, '"',
                ' "', file_mod, '"')
   
  write(cmd,file="script.txt",append=TRUE)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

Now’s the time to step out of R and into GDAL

1. In the start menu, search for OSGeo4W Shell and open it. You should see a command-line style interface.
2. In your working directory, find and open script.txt. It should have a bunch of commands in it.
3. Copy all the commands in script.txt and paste them right into the OSGeo4W Shell. Press enter to run the commands.
4. If you have many shapefiles, this might take a while. Wait until everything has finished running before moving to the next step.

Now, hop back into R.

Move the files we just created back into the working directory.

files <- list.files(tempdir(), pattern = "_reproj.tif")
for (i in 1:length(files)) {
  file.rename(from=paste0(tempdir(), "/", files[i]),to= paste0("./reproj/", files[i]))
}

1
2
3
4
5

Now get the average color of each image. Here I use a slightly different method than I did for the plots shown above. While writing up this code, I found a neat trick (resizing the image down to 1 px) that works faster than what I did before (averaging the R, G, and B values of each pixel).

files <- list.files("./reproj/")
colors <- NULL
 
for (i in 1:length(files)) {
  
  img <- image_read(paste0("./reproj/", files[i]))
  img <- image_transparent(img, "#ffffff", fuzz = 0)
 
  
  
  img <- image_scale(img,"1x1!")
   
  
  avgcolor <- image_data(img, channels = "rgba")[1:3] %>% paste0(collapse="")
   
  
  colors <- bind_rows(colors, data.frame(name = gsub("_reproj.tif", "", files[i]), color = paste0("#", avgcolor)))
}
 

write.csv(colors, "colors.csv", row.names = F)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

Et voilà! A lovely csv with the average color of each of your regions.

And if you’d like to plot it all pretty, here’s how I did it in R. For this example, I’m only going to plot the lower 48 states.

install.packages(c('sf', 'ggplot2', 'ggthemes', 'dplyr'))
 

library(sf)
library(ggplot2)
library(ggthemes)
library(dplyr)
 

setwd("C:/Users/Erin/Documents/DataViz/AvgColor/Temp/")
 

shapes <- read_sf("cb_2018_us_state_20m.shp")
colors <- read.csv("colors.csv")
 

shapes <- inner_join(shapes, colors, by = c("NAME" = "name")) 
 


shapes <- st_transform(shapes, 2163)
 

shapes <- subset(shapes, !(STUSPS %in% c("AK", "HI", "PR")))
 

ggplot() + theme_map() +
  geom_sf(data = shapes, aes(fill = I(color)), color = alpha("#fffbf2", .1)) +
  theme(plot.background = element_rect(fill = '#fffbf2', colour = '#fffbf2'))
 

ggsave("avgcolors.png", width = 7, height = 9, units = "in", dpi = 500)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

---

Erin Davis

Originally posted on the 26th of June 2021 at https://erdavis.com/ (opens new window)

Link to the article: https://erdavis.com/2021/06/26/average-colors-of-the-world/ (opens new window)

Sentinel-2 cloudless - https://s2maps.eu (opens new window) by EOX IT Services GmbH (Contains modified Copernicus Sentinel data 2019)

---

Media coverage and other references:

from Hacker News https://ift.tt/3eoezKt