Как правильно добавить карту в растровое изображение в R

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

Чтобы воспроизвести эту проблему, вот ссылка на папку с файлом, с которым я работаю: https://www.dropbox.com/s/e8pwgmnhvw4s0nf/sst.nc4?dl=0

Вот код, который я разработал до сих пор;

library(ncdf4)
library(raster)
library(mapdata)
library(mapproj)
library(rgeos)
library(ggplot2)

Через ncdf4, rasterToPoints, map_data и ggplot2

eight = nc_open("Downloads/sst.nc4")
sst = ncvar_get(eight, "sst")
sst = raster(sst)
sst = t(flip(sst, 1)) # have to orient the data properly

# extract the dimensions and set the extent
lat.min = min(eight$dim$lat$vals)
lat.max = max(eight$dim$lat$vals)
lon.min = min(eight$dim$lon$vals)
lon.max = max(eight$dim$lon$vals)
sst = setExtent(sst, ext = c(lon.min, lon.max, lat.min, lat.max))

# provide proper projection
crs(sst) = "+init=epsg:4326"

# convert raster to points
sst.p <- rasterToPoints(sst)
df <- data.frame(sst.p)
colnames(df) <- c("Longitude", "Latitude", "sst")
usa = map_data("usa")
ggplot(data=df, aes(y=Latitude, x=Longitude)) +
  geom_raster(aes(fill=sst)) +
  theme_bw() +
  coord_equal() +
  scale_fill_gradient("SST (Celsius)", limits=c(0,35)) +
  geom_polygon(data = usa, aes(x=long, y = lat, group = group)) + 

  theme(axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=16, angle=90),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right",
        legend.key = element_blank()
  )

ggplot.png

Через растр, данные maptools, SP Transform и базовое построение

#read in the data
sst = raster("Downloads/sst.nc4",  varname = "sst", stopIfNotEqualSpaced=FALSE)

# get world map data
data("wrld_simpl", package="maptools")

## Crop to the desired extent, then plot
newext <- c(lon.min, lon.max, lat.min, lat.max)
out <- crop(wrld_simpl, newext)

#transform to proper CRS
out = spTransform(out, "+init=epsg:4326")

#plot
plot(out, col="khaki", bg="azure2")
plot(sst, add = T)

baseplot.png

- Проекция, которую я использую для этих пространственных данных, EPSG:4326

-Вот фрагмент XML, определяющий выходную проекцию sst.nc4

<crs>PROJCS["Mercator_1SP / World Geodetic System 1984",
         GEOGCS["World Geodetic System 1984",
         DATUM["World Geodetic System 1984",
         SPHEROID["WGS 84", 6378135.0, 298.257223563, AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
         UNIT["degree", 0.017453292519943295],
         AXIS["Geodetic longitude", EAST],
         AXIS["Geodetic latitude", NORTH]],
         PROJECTION["Mercator_1SP"],
         PARAMETER["latitude_of_origin", 0.0],
         PARAMETER["central_meridian", 0.0],
         PARAMETER["scale_factor", 1.0],
         PARAMETER["false_easting", 0.0],
         PARAMETER["false_northing", 0.0],
         UNIT["m", 1.0],
         AXIS["Easting", EAST],
         AXIS["Northing", NORTH]]</crs>

Я также пытался использовать функцию map() с аргументом projection mapproj, но, похоже, у него нет проекции псевдомеркатора в качестве опции.


r geospatial projection satellite-image
person James    schedule 14.11.2017    source источник
comment
Можете ли вы получить это правильно в arcmap или Qgis? Мне кажется, у вас проблемы с проекцией. Вы получаете данные от AQUA/MODIS?   -  person M--    schedule 15.11.2017
comment
@Masoud Я не пробовал arcmap или Qgis. Знаете ли вы учебник для этих пакетов? Это данные AQUA/MODIS.   -  person James    schedule 15.11.2017
comment
Это не пакеты R. Это программное обеспечение gis. У НАСА есть программа под названием Panoply (или похожее название). Установите его и загрузите в него файл netcdf. Если это сработало, вам нужно изменить проекцию вашей базовой карты в R, чтобы получить желаемый результат.   -  person M--    schedule 15.11.2017

Ответы (1)


arrow_upward
2
arrow_downward

Это немного сбивает с толку. Обычно самым простым подходом было бы

sst = raster("sst.nc4",  varname = "sst")

но для этого файла это дает эту ошибку:

"cells are not equally spaced; you should extract values as points"

Итак, давайте сделаем это:

library(ncdf4)
library(raster)
library(maptools)

d <- nc_open("sst.nc4")
sst <- ncvar_get(d, "sst")
lon <- ncvar_get(d, "lon")
lat <- ncvar_get(d, "lat")
nc_close(d)

xy <- cbind(rep(lon, length(lat)), rep(lat, each=length(lon)))

Объедините и удалите значения NA (около половины ячеек...

xyv <- na.omit(cbind(xy, as.vector(sst)))

Настройте RasterLayer с разрешением, достаточным для ваших целей, и растеризуйте точки

 r <- raster(extent(range(lon), range(lat)), res=1/6)
 r <- rasterize(xyv[, 1:2], r, xyv[,3], fun=mean)

Сюжет

data(wrld_simpl)
w <- crop(wrld_simpl, r)

plot(r)
plot(w, col='gray', add=TRUE)
person Robert Hijmans    schedule 14.11.2017
comment
это работает! Будет ли этот метод работать без исключения NA? Я хотел бы оставить NA, если это возможно. Мои первоначальные попытки отредактировать это не увенчались успехом - person James; 15.11.2017
comment
Вы можете оставить NA, но это не изменит результатов. - person Robert Hijmans; 15.11.2017
comment
это из-за растеризации, которая, кажется, имеет какую-то среднюю функциональность для заполнения пробелов? - person James; 15.11.2017