Репроекция океанографических данных в растр?
Относительно моего предыдущего вопроса Чтение / открытие данных Oceancolor из MODIS в формате nc. Я попытался прочитать загрузку данных ncdf oceancolor здесь, используя ncdf4 package
с этим кодом
library(ncdf4)
nc <- nc_open('A2014325053500.L2_LAC_SST.nc')
# Get data for each variabel name as a large matrix
sst <- ncvar_get(nc, varid = 'geophysical_data/sst')
lon <- ncvar_get(nc, varid = 'navigation_data/longitude')
lat <- ncvar_get(nc, varid = 'navigation_data/latitude')
чем эти данные, я конвертирую в растр. Но мне не удалось перепроектировать растр figure 1
, Изображение должно быть как figure 2
library(sp)
library(raster)
r <- raster(sst,
xmn=min(lat), xmx=max(lat),
ymn=min(lon), ymx=max(lon),
crs=CRS('+proj=longlat +ellps=WGS84 +datum=WGS84'))
proj='+proj=longlat +ellps=WGS84 +datum=WGS84'
r2 <- projectRaster(r, crs=proj, method = 'ngb')
plot(r2)
Исходя из первого предположения jbaums, результат неправильно корректируется figure 3
Я просто понимаю, что что-то не так с результатом, когда я увеличиваю разрешение. В западной и восточной части становятся белыми (нет данных). Я использовал этот код ниже:
r <- rasterize(df, raster(extent(df), res=0.0103011), 'sst', fun=mean)
1 ответ
В их текущей проекции координаты не попадают в регулярную сетку.
Одним из способов решения этой проблемы является представление координат в виде SpatialPoints
объект, а затем, если вам нужен растр, rasterize
Это.
library(ncdf4)
library(raster)
# Read in the data
nc <- nc_open('~/../Downloads/A2014325053500.L2_LAC_SST.nc')
sst <- ncvar_get(nc, varid = 'geophysical_data/sst')
lon <- ncvar_get(nc, varid = 'navigation_data/longitude')
lat <- ncvar_get(nc, varid = 'navigation_data/latitude')
nc_close(nc)
# Create a SpatialPointsDataFrame
p <- data.frame(lon=c(lon), lat=c(lat), sst=c(sst))
coordinates(p) <- ~lon+lat
proj4string(p) <- '+init=epsg:4326'
# Rasterize with appropriate resolution and aggregation function
r <- rasterize(p, raster(extent(p) * 1.04, res=0.05), 'sst', fun=mean)
# for high resolution look at ?gdalUtils::gdal_rasterize for efficiency
Вот сюжет:
library(rasterVis)
levelplot(r, at=seq(-75, 75, length=100), margin=FALSE, par.settings=BuRdTheme,
colorkey=list(height=0.6), main='MODIS Sea Surface Temperature')