why is tmap function within for loop that also has left join not working

Question

I have a shapefile "grid" that is basically a grid map of US and Canada. Each grid has its own specific ID (grid_code). I also have a csv file that has the relative abundance values of bird species that correlate to respective grids. Because each species has their unique distribution, the species table only includes a subset of all the available grid codes.

I am trying to create a thematic map showing the relative abundance of each species. For this, I am first joining the csv file to the attribute table of the grid shapefile using leftjoin function, where grid_code is the common column. After specifying the color palette for the tmap, I am creating a png file for each species' relative abundance.

The following code works fine when I manually put in a value for i; however, when I run the entire code, it stops after the leftjoin function, and creates a blank png file for each species. I have tried several things but am having no luck fixing this. Is it just a matter of placement of codes or is there something majorly wrong with my for-loop?

Thanks.

# reading in csv file with relative abundance info
rel_ab <- read.csv("ra_test.csv", h=T, stringsAsFactors = F)

# specifying color palette for tmap
colfunc <- colorRampPalette(c("lightpink", "red"))

# running for loop code where I'm reading and reprojecting the grid
# shapefile each time, then joining it to the species table, and finally 
# plotting the tmap

# only want to work with one species at a time, so first:
ra <- unique(rel_ab$species)

for (i in ra){
   species <- rel_ab[which(rel_ab$species == i),]
   grid <- readOGR(dsn = "grid_clip", "grid") # reading grid shapefile
   grid83<- spTransform(grid, CRS("+init=epsg:5070")) # transorming shapefile
   grid83@data <- left_join(grid83@data, species) 

# grid83 now contains a column that includes all the grid codes PLUS a 
# column for relative abundance values (V2) for species i. Again, because a
# species is not found in every part of US/CANADA, there are a lot of V2
# rows that are NA.

# Here is what grid83@data looks like:

head(grid83@data)

AREA PERIMETER TEMP_ TEMP_ID GRID_CODE BCR STATE STATE_CODE STRATA_ KMSQ  X V2 V3 species
1 463738112  87471.12  6548   92624      6981   5   ALK          3     S94  461 NA NA NA    <NA>
2 463749632  87485.10  6551   92625      6982   5   ALK          3     S94  461 NA NA NA    <NA>
3 463766720  87499.32  6554   92626      6983   4    YT         93     S68  461 NA NA NA    <NA>

   png(width = 1000, height = 1000, filename = paste("RA_maps/",i,".png", sep=""))

   tm_shape(grid83)+
    tm_fill("V2", title = "BBS Summer Distribution", textNA = "None counted", 
  style = "fixed", breaks = c(Inf, 100, 30, 10, 3, 1, 0.05), colorNA = "white", palette = colfunc(6))+
    tm_shape(uscan83)+ 
    tm_borders(lwd = 1.5)+
   tm_shape(gray83)+
    tm_fill(col="grey85")+
    tm_borders()+
   tm_layout(
     inner.margins = c(.02, .02, .02, .02),
     legend.title.size = 1.8,
     legend.text.size = 1,
     legend.bg.color = "white",
     legend.bg.alpha = 0,
     legend.width = 0.22)

   dev.off()}

could you paste some of your data please? dput(rel_ab) or if thats too big then dput(head(rel_ab)) in your post? — andrnev, Jan 29 '16 at 23:44
@andrnev, grid83 is the shapefile and grid83@data outputs the attribute table of the shapefile. grid83@data after the leftjoin looks like this: [removed] because of formatting issues, I'm going to add the output on the original question. — wallflower, Jan 30 '16 at 00:02
Not too sure whether this will work but could you modify the following grid83@data <- asS4(left_join(asS3(grid83@data), species)) and check? I am not too sure whether the left_join function works on S4 class objects — andrnev, Jan 30 '16 at 00:37
@andrnev - hmm, that didn't work. Got this error msg: Error in asS3(grid83@data, species) : (list) object cannot be coerced to type 'logical' — wallflower, Jan 30 '16 at 00:47
@andrnev - I just double checked and there are no NAs in ra. Also, the code runs fine when I put in an "i" value myself. I get a nice tmap and its saved in the right folder and everything. But as soon as I run the code as a for loop it stops right after the leftjoin. It does create a png file in the right folder, but its empty. That's why I'm wondering whether its just an issue with placement of codes rather than the code itself. — wallflower, Jan 30 '16 at 01:10
Yep. Could you try to traceback the problem. Have a look at this for debugging steps, http://www.biostat.jhsph.edu/~hji/courses/statcomputing/Debug.pdf — andrnev, Jan 30 '16 at 09:32

score 3 · Answer 1 · answered Feb 08 '16 at 13:14

You'll have to print the tmap plot explicitly if you run it inside a for loop. Alternatively, there is a save_tmap in the package that could be useful. See the following code chunk, which also contains some other tmap function. Since I don't have the shape and the data, I cannot test it myself.

for (i in ra){
   species <- rel_ab[which(rel_ab$species == i),]
   grid83 <- read_shape("grid_clip/grid.shp", projection = 5070)
   grid83 <- append_data(grid83, species) 

   m <- tm_shape(grid83)+
    tm_fill("V2", title = "BBS Summer Distribution", textNA = "None counted", 
  style = "fixed", breaks = c(Inf, 100, 30, 10, 3, 1, 0.05), colorNA = "white", palette = colfunc(6))+
    tm_shape(uscan83)+ 
    tm_borders(lwd = 1.5)+
   tm_shape(gray83)+
    tm_fill(col="grey85")+
    tm_borders()+
   tm_layout(
     inner.margins = c(.02, .02, .02, .02),
     legend.title.size = 1.8,
     legend.text.size = 1,
     legend.bg.color = "white",
     legend.bg.alpha = 0,
     legend.width = 0.22)

   save_tmap(m, width = 1000, height = 1000, units="px", filename = paste("RA_maps/",i,".png", sep=""))
}

why is tmap function within for loop that also has left join not working

1 Answers1

Linked