почему функция tmap внутри цикла for, которая также оставила соединение, не работает

У меня есть "сетка" из шейп-файла, которая в основном представляет собой сеточную карту США и Канады. Каждая сетка имеет свой собственный идентификатор (grid_code). У меня также есть файл csv, в котором есть значения относительной численности видов птиц, которые коррелируют с соответствующими сетками. Поскольку каждый вид имеет свое уникальное распространение, таблица видов включает только подмножество всех доступных кодов сетки.

Я пытаюсь создать тематическую карту, показывающую относительную численность каждого вида. Для этого я сначала присоединяю файл csv к таблице атрибутов шейп-файла сетки с помощью функции leftjoin, где grid_code является общим столбцом. После определения цветовой палитры для tmap я создаю файл png для относительной численности каждого вида.

Следующий код отлично работает, когда я вручную ввожу значение для i; однако, когда я запускаю весь код, он останавливается после функции leftjoin и создает пустой файл png для каждого вида. Я пробовал несколько вещей, но мне не удалось это исправить. Это просто вопрос размещения кодов или что-то не так с моим циклом for?

Спасибо.

# 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()}

r left-join png gis tmap
person wallflower    schedule 29.01.2016    source источник
comment
не могли бы вы вставить некоторые из ваших данных? dput (rel_ab) или, если он слишком большой, то dput (head (rel_ab)) в вашем сообщении?   -  person andrnev    schedule 30.01.2016
comment
И что такое grid83 @ data?   -  person andrnev    schedule 30.01.2016
comment
@andrnev, grid83 - шейп-файл, а grid83 @ data - выводит таблицу атрибутов шейп-файла. grid83 @ data после leftjoin выглядит так: [удалено] из-за проблем с форматированием, я собираюсь добавить вывод по исходному вопросу.   -  person wallflower    schedule 30.01.2016
comment
Не уверен, что это сработает, но не могли бы вы изменить следующую grid83 @ data ‹- asS4 (left_join (asS3 (grid83 @ data), sizes)) и проверить? Я не слишком уверен, работает ли функция left_join с объектами класса S4.   -  person andrnev    schedule 30.01.2016
comment
@andrnev - хм, не сработало. Получил эту ошибку msg: Ошибка в asS3 (grid83 @ data, sizes): (list) объект не может быть принудительно приведен к типу 'logic'   -  person wallflower    schedule 30.01.2016
comment
Может быть, некоторые ценности в ra принадлежат НС?   -  person andrnev    schedule 30.01.2016
comment
@andrnev - Я только что перепроверил, и в ra нет НАС. Кроме того, код работает нормально, когда я сам ввожу значение i. Я получаю хороший tmap, сохраняю его в нужной папке и все такое. Но как только я запускаю код как цикл for, он останавливается сразу после leftjoin. Он создает файл png в правой папке, но он пустой. Вот почему мне интересно, не проблема ли это просто в размещении кодов, а не в самом коде.   -  person wallflower    schedule 30.01.2016
comment
Ага. Не могли бы вы попытаться отследить проблему. Взгляните на это для шагов отладки, biostat.jhsph.edu/ ~ hji / курсы / statcomputing / Debug.pdf   -  person andrnev    schedule 30.01.2016

Ответы (1)


arrow_upward
2
arrow_downward

Вам придется print график tmap явно, если вы запустите его внутри цикла for. В качестве альтернативы в пакете есть save_tmap, который может быть полезен. См. Следующий фрагмент кода, который также содержит некоторые другие функции tmap. Поскольку у меня нет формы и данных, я не могу проверить это сам.

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=""))
}
person Martijn Tennekes    schedule 08.02.2016