Переменная из файла netcdf выходит перевернутой

Я скачал файл nc из

f=open.ncdf("file.nc")
[1] "file Lfile.nc has  2 dimensions:"
[1] "Longitude   Size: 1440"
[1] "Latitude   Size: 720"
[1] "------------------------"
[1] "file filr.nc has   8 variables:"
[1] "short ts[Latitude,Longitude]  Longname:Skin Temperature (2mm) Missval:NA"

Затем я хотел работать с переменной ground_moisture_c

A = get.var.ncdf(nc=f,varid="soil_moisture_c",verbose=TRUE)

Я тогда заговор с image(A), Я получил карту, показанную ниже, я даже перенес ее image(t(a)) но это было изменено в другом направлении, а не так, как должно быть. В любом случае, чтобы узнать, что не так, я использовал средство просмотра Netcdf Panoply, и карта была правильно построена, как вы можете видеть ниже.

Источник

3 ответа

Решение

В дополнение к гораздо лучшему решению @mdsumner, здесь есть способ сделать это с помощью библиотеки. ncdf только.

library(ncdf)
f <- open.ncdf("LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040101.nc")
A <- get.var.ncdf(nf,"soil_moisture_c")

Все, что вам нужно, это найти ваши размеры, чтобы иметь согласованную ось X и Y. Если вы посмотрите на размеры ваших объектов netCDF, то вот что вы видите:

str(f$dim)
List of 2
 $ Longitude:List of 8
  ..$ name         : chr "Longitude"
  ..$ len          : int 1440
  ..$ unlim        : logi FALSE
  ..$ id           : int 1
  ..$ dimvarid     : num 2
  ..$ units        : chr "degrees_east"
  ..$ vals         : num [1:1440(1d)] -180 -180 -179 -179 -179 ...
  ..$ create_dimvar: logi TRUE
  ..- attr(*, "class")= chr "dim.ncdf"
 $ Latitude :List of 8
  ..$ name         : chr "Latitude"
  ..$ len          : int 720
  ..$ unlim        : logi FALSE
  ..$ id           : int 2
  ..$ dimvarid     : num 1
  ..$ units        : chr "degrees_north"
  ..$ vals         : num [1:720(1d)] 89.9 89.6 89.4 89.1 88.9 ...
  ..$ create_dimvar: logi TRUE
  ..- attr(*, "class")= chr "dim.ncdf"

Следовательно, ваши размеры:

 f$dim$Longitude$vals -> Longitude
 f$dim$Latitude$vals -> Latitude

Теперь ваш Latitude идет от 90 до -90 вместо противоположного, которое image Предпочитает, так:

 Latitude <- rev(Latitude)
 A <- A[nrow(A):1,]

Наконец, как вы заметили, x и y вашего объекта A перевернуты, поэтому вам нужно переставить его, и ваши значения NA по некоторым причинам представлены по значению -32767:

A[A==-32767] <- NA
A <- t(A)

И наконец сюжет:

image(Longitude, Latitude, A)
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)

Изменить: чтобы сделать это на ваших 31 файлов, давайте назовем ваш вектор имен файлов ncfiles а также yourpath каталог, где вы их сохранили (для простоты я собираюсь предположить, что ваша переменная всегда вызывается soil_moisture_c и ваши НС всегда -32767):

ncfiles
 [1] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040101.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040102.nc"
 [3] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040103.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040104.nc"
 [5] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040105.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040106.nc"
 [7] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040107.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040108.nc"
 [9] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040109.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040110.nc"
[11] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040111.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040112.nc"
[13] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040113.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040114.nc"
[15] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040115.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040116.nc"
[17] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040117.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040118.nc"
[19] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040119.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040120.nc"
[21] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040121.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040122.nc"
[23] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040123.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040124.nc"
[25] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040125.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040126.nc"
[27] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040127.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040128.nc"
[29] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040129.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040130.nc"
[31] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040131.nc"

yourpath
 [1] "C:\\Users"

library(ncdf)
library(maptools)
data(wrld_simpl)
for(i in 1:length(ncfiles)){
    f <- open.ncdf(paste(yourpath,ncfiles[i], sep="\\"))
    A <- get.var.ncdf(f,"soil_moisture_c")
    f$dim$Longitude$vals -> Longitude
    f$dim$Latitude$vals -> Latitude
    Latitude <- rev(Latitude)
    A <- A[nrow(A):1,]
    A[A==-32767] <- NA
    A <- t(A)
    close.ncdf(f) # this is the important part
    png(paste(gsub("\\.nc","",ncfiles[i]), "\\.png", sep="")) # or any other device such as pdf, jpg...
    image(Longitude, Latitude, A)
    plot(wrld_simpl, add = TRUE)
    dev.off()
    }
Источник

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

library(raster) ## requires ncdf package for this file  
d <- raster("LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T185959Z_20040114.nc", varname = "soil_moisture_c")

(Я использовал другой файл для вашего, но он должен работать так же).

Оказывается, что даже растр не получает это право без работы, но это действительно облегчает исправление:

d <-  flip(t(d), direction = "x")

Это транспонировало данные и переключалось вокруг "x" (долгота), сохраняя географическую привязку в исходном контексте.

Подготовьте карту с помощью maptools, чтобы проверить:

plot(d)

library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)

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

сюжет изображения из растра

Источник

Вы также можете просто инвертировать широты из командной строки, используя CDO:

cdo invertlat file.nc file_inverted.nc
Источник
Другие вопросы по тегам