Извлечение матрицы из NetCDF и преобразование ее в растр - проблемы со строками - R - PullRequest
Новая
       60

Извлечение матрицы из NetCDF и преобразование ее в растр - проблемы со строками - R

0 голосов
/ 29 апреля 2020

У меня есть следующая матрица, полученная из файла NetCDF, который я пытаюсь преобразовать в растр. Я знаю, что матрица имеет проекцию WGS84.

Я обнаружил одну проблему, исправленную в приведенном ниже коде - интервал между широтами lat в NetCDF и, следовательно, в производной матрице clim_ncdf не был равен. Тем не менее, растр EI_adj не проецируется точно, как это должно быть после преобразования из матрицы, и все пространства смещены на юг. Эта проблема сводит меня с ума - у кого-нибудь есть идеи, как это исправить? Исходные файлы (NetCDF и границы мировых администраторов) можно загрузить с здесь . 

library(raster)
library(ncdf4)
library(lattice)

# Choose variable name
dname <- c("GI")

clim_ncdf <- nc_open("NetCDF_GI.nc")

lon <- ncvar_get(clim_ncdf,"Longitude")
head(lon)
lat <- ncvar_get(clim_ncdf,"Latitude")

# Latitudes have spacing of 0.5 except in two instances:
lat[55:60]
lat[58:59] # problematic ones
# Create a new latitude vector with equal spacing for corrected matrix
nlat <- seq(min(lat),max(lat),0.5)

EI1 <- ncvar_get(clim_ncdf,dname[1])
# This needs to be rotated
rotate <- function(x) t(apply(x, 2, rev))
EI <- rotate(rotate(rotate(EI1)))

# Now adjust EI for the problematic lats:
rows_m_reps <- rep(1,nrow(EI))
rows_m_reps[58] <- 2
rows_m_reps[59] <- 10

# Replicating corresponding rows so we can now have equal latitude distancing
EI_adj <- EI[rep(1:nrow(EI), rows_m_reps), ] 

EIr_adj <- raster(EI_adj,xmn=min(lon), xmx=max(lon),ymn=min(nlat),ymx=max(nlat),
                  crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(EIr_adj)

# Add admin layer
library(rgdal)
world_admin <- readOGR("Countries_WGS84.shp")
plot(world_admin, add = TRUE)

enter image description here

1 Ответ

2 голосов
/ 29 апреля 2020

Ваши данные 

library(ncdf4)
library(raster)
library(maptools)
data(wrld_simpl)

clim_ncdf <- nc_open("NetCDF_GI.nc")
lon <- ncvar_get(clim_ncdf,"Longitude")
lat <- ncvar_get(clim_ncdf,"Latitude")
v <- ncvar_get(clim_ncdf, "GI")

Это все широты. 

lat2 <- seq(min(lat),max(lat),0.5)

Создайте RasterLayer и соответствующую матрицу с NA s 

e <- extent(min(lon)-0.25, max(lon)+0.25, min(lat)-0.25, max(lat)+0.25)
r <- raster(nrow=length(lat2), ncol=length(lon), ext=e)
m <- matrix(NA, nrow=length(lat2), ncol=length(lon))

Теперь поверните и присвойте значения правильным строкам для матрицы m

vv <- t(v[,ncol(v):1])
i <- rowFromY(r, rev(as.vector(lat)))
m[i,] <- vv

И присвоить m растровому слою 

values(r) <- m
image(r)
lines(wrld_simpl)
...