отаю со спутниковыми данными, организованными на нерегулярной двумерной сетке, размеры которой - линия развертки (по размеру дорожки) и пиксель земли (по измерению дорожки). Информация о широте и долготе для каждого центрального пикселя хранится во вспомогательных переменных координат, а также в парах координат четырех углов (координаты широты и долготы указаны в эталонном эллипсоиде WGS84). Данные хранятся в файлах netCDF4.

Я пытаюсь эффективно нанести эти файлы (и, возможно, комбинацию файлов - следующий шаг!) На проецируемую карту.

Мой подход до сих пор вдохновленДжереми Войси ответ на этовопросбыло создание фрейма данных, который связывает интересующую меня переменную с границами пикселей, и использоватьggplot2 с участиемgeom_polygon для фактического сюжета.

Позвольте мне проиллюстрировать мой рабочий процесс и заранее извиняюсь за наивный подход: я только начал программировать на R с недели или двух.

Запись

Чтобы полностью воспроизвести проблему:
1. скачать два кадра данных:so2df.Rda (22M) иpixel_corners.Rda (26М)
2. загрузить их в вашу среду, например,

so2df <- readRDS(file="so2df.Rda")
pixel_corners <- readRDS(file="pixel_corners.Rda")

перейдите к шагу «Объединить кадры данных».Начальная настройка

Я собираюсь прочитать данные и границы широты / долготы из моего файла.

library(ncdf4)
library(ggplot2)
library(ggmap) 
# set path and filename
ncpath <- "/Users/stefano/src/s5p/products/e1dataset/L2__SO2/"
ncname <- "S5P_OFFL_L2__SO2____20171128T234133_20171129T003956_00661_01_022943_00000000T000000"  
ncfname <- paste(ncpath, ncname, ".nc", sep="")
nc <- nc_open(ncfname)

# save fill value and multiplication factors
mfactor = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column", 
                    "multiplication_factor_to_convert_to_DU")
fillvalue = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column", 
                      "_FillValue")

# read the SO2 total column variable
so2tc <- ncvar_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column")

# read lat/lon of centre pixels
lat <- ncvar_get(nc, "PRODUCT/latitude")
lon <- ncvar_get(nc, "PRODUCT/longitude")

# read latitude and longitude bounds
lat_bounds <- ncvar_get(nc, "GEOLOCATIONS/latitude_bounds")
lon_bounds <- ncvar_get(nc, "GEOLOCATIONS/longitude_bounds")

# close the file
nc_close(nc)
dim(so2tc)
## [1]  450 3244

Таким образом, для этого файла / прохода есть 450 наземных пикселей для каждой из 3244 строк развертки.

Создайте кадры данных

Здесь я создаю два кадра данных, один для значений, с некоторой последующей обработкой, а другой для границ широты / долготы, а затем объединяю два кадра данных.

so2df <- data.frame(lat=as.vector(lat), lon=as.vector(lon), so2tc=as.vector(so2tc))
# add id for each pixel
so2df$id <- row.names(so2df)
# convert to DU
so2df$so2tc <- so2df$so2tc*as.numeric(mfactor$value)
# replace fill values with NA
so2df$so2tc[so2df$so2tc == fillvalue] <- NA
saveRDS(so2df, file="so2df.Rda")
summary(so2df)

##       lat              lon              so2tc              id           
##  Min.   :-89.97   Min.   :-180.00   Min.   :-821.33   Length:1459800    
##  1st Qu.:-62.29   1st Qu.:-163.30   1st Qu.:  -0.48   Class :character  
##  Median :-19.86   Median :-150.46   Median :  -0.08   Mode  :character  
##  Mean   :-13.87   Mean   : -90.72   Mean   :  -1.43                     
##  3rd Qu.: 31.26   3rd Qu.: -27.06   3rd Qu.:   0.26                     
##  Max.   : 83.37   Max.   : 180.00   Max.   :3015.55                     
##                                     NA's   :200864

Я сохранил этот фрейм данных какso2df.Rda Вот (22М).

num_points = dim(lat_bounds)[1]
pixel_corners <- data.frame(lat_bounds=as.vector(lat_bounds), lon_bounds=as.vector(lon_bounds))
# create id column by replicating pixel's id for each of the 4 corner points
pixel_corners$id <- rep(so2df$id, each=num_points)
saveRDS(pixel_corners, file="pixel_corners.Rda")
summary(pixel_corners)


##    lat_bounds       lon_bounds           id           
##  Min.   :-89.96   Min.   :-180.00   Length:5839200    
##  1st Qu.:-62.29   1st Qu.:-163.30   Class :character  
##  Median :-19.86   Median :-150.46   Mode  :character  
##  Mean   :-13.87   Mean   : -90.72                     
##  3rd Qu.: 31.26   3rd Qu.: -27.06                     
##  Max.   : 83.40   Max.   : 180.00

Как и ожидалось, фрейм данных границ широты и долготы в четыре раза больше фрейма данных значения (четыре точки для каждого пикселя / значения).
Я сохранил этот фрейм данных какpixel_corners.Rda Вот (26М).

Объединить кадры данных

Затем я объединяю два фрейма данных по id:

start_time <- Sys.time()
so2df <- merge(pixel_corners, so2df, by=c("id"))
time_taken <- Sys.time() - start_time
print(paste(dim(so2df)[1], "rows merged in", time_taken, "seconds"))

## [1] "5839200 rows merged in 42.4763631820679 seconds"

Как видите, это процесс, который сильно загружает процессор. Интересно, что произойдет, если я буду работать с 15 файлами одновременно (глобальный охват).

Построение данных

Теперь, когда у меня есть углы пикселей, связанные со значением пикселя, я могу легко построить их. Обычно меня интересует конкретная область орбиты, поэтому я сделал функцию, которая подставляет входной фрейм данных перед его построением:

PlotRegion <- function(so2df, latlon, title) {
  # Plot the given dataset over a geographic region.
  #
  # Args:
  #   df: The dataset, should include the no2tc, lat, lon columns
  #   latlon: A vector of four values identifying the botton-left and top-right corners 
  #           c(latmin, latmax, lonmin, lonmax)
  #   title: The plot title

  # subset the data frame first
  df_sub <- subset(so2df, lat>latlon[1] & lat<latlon[2] & lon>latlon[3] & lon<latlon[4])

  subtitle = paste("#Pixel =", dim(df_sub)[1], "- Data min =", 
                   formatC(min(df_sub$so2tc, na.rm=T), format="e", digits=2), "max =", 
                   formatC(max(df_sub$so2tc, na.rm=T), format="e", digits=2))

  ggplot(df_sub) + 
    geom_polygon(aes(y=lat_bounds, x=lon_bounds, fill=so2tc, group=id), alpha=0.8) +
    borders('world', xlim=range(df_sub$lon), ylim=range(df_sub$lat), 
            colour='gray20', size=.2) + 
    theme_light() + 
    theme(panel.ontop=TRUE, panel.background=element_blank()) +
    scale_fill_distiller(palette='Spectral') +
    coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
    labs(title=title, subtitle=subtitle, 
         x="Longitude", y="Latitude", 
         fill=expression(DU)) 
}

Затем я вызываю свою функцию в интересующей области, например, давайте посмотрим, что происходит на Гавайях:

latlon = c(17.5, 22.5, -160, -154)
PlotRegion(so2df, latlon, expression(SO[2]~total~vertical~column))

Вот они, мои пиксели, и, похоже, шлейф SO2 от Мауна Лоа. Пожалуйста, пока игнорируйте отрицательные значения. Как вы можете видеть, область пикселей изменяется к краю полосы (другая схема биннинга).

Я попытался показать тот же график на картах Google, используя ggmap:

PlotRegionMap <- function(so2df, latlon, title) {
  # Plot the given dataset over a geographic region.
  #
  # Args:
  #   df: The dataset, should include the no2tc, lat, lon columns
  #   latlon: A vector of four values identifying the botton-left and top-right corners 
  #           c(latmin, latmax, lonmin, lonmax)
  #   title: The plot title

  # subset the data frame first
  df_sub <- subset(so2df, lat>latlon[1] & lat<latlon[2] & lon>latlon[3] & lon<latlon[4])

  subtitle = paste("#Pixel =", dim(df_sub)[1], "Data min =", formatC(min(df_sub$so2tc, na.rm=T), format="e", digits=2), 
                   "max =", formatC(max(df_sub$so2tc, na.rm=T), format="e", digits=2))
  base_map <- get_map(location = c(lon = (latlon[4]+latlon[3])/2, lat = (latlon[1]+latlon[2])/2), zoom = 7, maptype="terrain", color="bw")

  ggmap(base_map, extent = "normal")  +
    geom_polygon(data=df_sub, aes(y=lat_bounds, x=lon_bounds,fill=so2tc, group=id),  alpha=0.5) +
    theme_light() + 
    theme(panel.ontop=TRUE, panel.background=element_blank()) +
    scale_fill_distiller(palette='Spectral') +
    coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
    labs(title=title, subtitle=subtitle, 
         x="Longitude", y="Latitude", 
         fill=expression(DU)) 

}

И вот что я получаю:

latlon = c(17.5, 22.5, -160, -154)
PlotRegionMap(so2df, latlon, expression(SO[2]~total~vertical~column))

ВопросовЕсть ли более эффективный способ решения этой проблемы? Я читал оsf пакет, и мне было интересно, смогу ли я определить датафрейм точек (значения + координаты центрального пикселя), и иметьsf автоматически выводить границы пикселей. Это избавило бы меня от необходимости полагаться на границы широты и долготы, определенные в моем исходном наборе данных, и на необходимость объединения их с моими значениями. Я мог бы принять потерю точности на переходных участках к краю полосы, в остальном сетка в значительной степени регулярна, каждый пиксель имеет размер 3.5x7 км ^ 2.Повысит ли производительность повторную привязку моих данных к регулярной сетке (как?), Возможно, путем объединения соседних пикселей? Я думал об использованииraster Пакет, который, как я понял, требует данных на обычной сетке. Это должно быть полезно в глобальном масштабе (например, графики по всей Европе), где мне не нужно отображать отдельные пиксели - на самом деле, я бы даже не смог их увидеть.Нужно ли перепроектировать мои данные при построении графика на карте Google?

[бонусные косметические вопросы]

Есть ли более элегантный способ подмоделировать мой фрейм данных для региона, обозначенного четырьмя угловыми точками?Как я могу изменить цветовую шкалу, чтобы более высокие значения выделялись относительно более низких значений? Я имел дело с масштабированием журнала с плохими результатами.