Оценка площади ячейки сетки по спутниковым данным с использованием стереографической полярной проекции

Я пытаюсь работать со спутниковыми данными из полярных регионов. Их можно скачать в формате .nc (netcdf, см. ниже). Они находятся в полярной стереографической проекции в регулярной сетке (https://nsidc.org/data/polar-stereo/ps_grids.html). Я хотел бы оценить площадь каждой ячейки, чтобы рассчитать площадь ледяного покрова, умножив долю ледяного покрова в каждой ячейке на площадь ячейки.

Я могу извлечь площадь ячейки, используя область из {растра}, только когда равномерная сетка находится в географических координатах (широта-долгота). Однако я не смог найти какую-либо функцию R, чтобы сделать это, когда сетка однородна в стереографических координатах.

Я новичок в пространственных данных и проекциях, и, возможно, я упускаю что-то важное.
Ниже некоторая важная информация и фрагмент кода с двумя испытаниями.

СООТВЕТСТВУЮЩИЕ ДАННЫЕ И ИСТОЧНИКИ ИНФОРМАЦИИ

Polar Watch https://polarwatch.noaa.gov/erddap/

Ссылка на данные о концентрации морского льда для загрузки данных

https://polarwatch.noaa.gov/erddap/griddap/nsidcCDRiceSQnhmday.nc?seaice_conc_monthly_cdr%5B(2019-12-16T00:00:00Z):1:(2019-12-16T00:00:00Z)%5D%5B(5837500.0):1:(-5337500.0)%5D%5B(-3837500.0):1:(3737500.0)%5D

Можно скачать в R с помощью

url1<- "https://polarwatch.noaa.gov/erddap/griddap/nsidcCDRiceSQnhmday.nc?seaice_conc_monthly_cdr[(2019-12-16T00:00:00Z):1:(2019-12-16T00:00:00Z)][(5837500.0):1:(-5337500.0)][(-3837500.0):1:(3737500.0)]"

download.file(url1, destfile='nsidcCDRiceSQnhmday_935c_47bd_a147.nc')

Данные для сетки Концентрация морского льда Сетка широты и долготы, запись климатических данных NOAA/NSIDC V3, Антарктика, 25 км ?широта%5B(4337500.0):1:(-3937500.0)%5D%5B(-3937500.0):1:(3937500.0)%5D,долгота%5B(4337500.0):1:(-3937500.0)%5D%5B(- 3937500.0):1:(3937500.0)%5D" rel="nofollow noreferrer">https://polarwatch.noaa.gov/erddap/griddap/nsidcCDRice_sh_grid.nc?latitude[(4337500.0):1:(-3937500.0)][ (-3937500.0):1:(3937500.0)],долгота[(4337500.0):1:(-3937500.0)][(-3937500.0):1:(3937500.0)]

url2<- "https://polarwatch.noaa.gov/erddap/griddap/nsidcCDRice_sh_grid.nc?latitude[(4337500.0):1:(-3937500.0)][(-3937500.0):1:(3937500.0)],longitude[(4337500.0):1:(-3937500.0)][(-3937500.0):1:(3937500.0)]"

download.file(url2, destfile="nsidcCDRice_sh_grid_c513_9d75_76c1.nc")

https://nsidc.org/data/polar-stereo/ps_grids.html Таблица 4. Проекция Южного полушария на основе WGS 1984 г. Широта истинного происхождения -70

НЕКОТОРЫЕ ИСПЫТАНИЯ

require(raster)
require(ncdf4)

br<-brick("nsidcCDRiceSQnhmday_935c_47bd_a147.nc")
projection(br)<- CRS("+init=epsg:3976 +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") #https://polarwatch.noaa.gov/tools-training/code-gallery/
br

res(br)

area(br)

который дает

Warning message:
In .local(x, ...) :
  This function is only useful for Raster* objects with a longitude/latitude coordinates

Использование сетки, предоставленной в источниках данных

IceFgrid<-nc_open("nsidcCDRice_sh_grid_c513_9d75_76c1.nc")

ygridLatLon <- ncvar_get(IceFgrid, varid="ygrid")
xgridLatLon <- ncvar_get(IceFgrid, varid="xgrid")
longitude <- ncvar_get(IceFgrid, varid="longitude")
latitude <- ncvar_get(IceFgrid, varid="latitude")
nc_close(IceFgrid)

dim(longitude) # Matrix with Longitude of each grid point.
length(xgridLatLon)
length(ygridLatLon)

## Try to convert to data frame and the to raster.

dims <- dim(longitude)
icemap.df <- data.frame(Longitude=array(longitude,dims[1]*dims[2]),
                        Latitude=array(latitude,dims[1]*dims[2]))
icemap.df$Seaice <- array(1,dims[1]*dims[2])# Here I use 1 for ice cover just

rast<-rasterFromXYZ(icemap.df, crs="+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs")

Но дает

Error in rasterFromXYZ(icemap.df, crs = "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs") : 
  x cell sizes are not regular

Любая помощь будет оценена.

Заранее спасибо Ангел


r area polar-coordinates satellite
person Angel Segura    schedule 06.05.2021    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Здесь я использую terra вместо raster

url1 <- "https://polarwatch.noaa.gov/erddap/griddap/nsidcCDRiceSQnhmday.nc?seaice_conc_monthly_cdr[(2019-12-16T00:00:00Z):1:(2019-12-16T00:00:00Z)][(5837500.0):1:(-5337500.0)][(-3837500.0):1:(3737500.0)]"

f <- 'nsidcCDRiceSQnhmday_935c_47bd_a147.nc'
download.file(url1, destfile='nsidcCDRiceSQnhmday_935c_47bd_a147.nc', mode="wb")
library(terra)
#terra version 1.3.4
r <- rast(f)
crs(r) <- "epsg:3976"

r
class       : SpatRaster 
dimensions  : 448, 304, 1  (nrow, ncol, nlyr)
resolution  : 25000, 25000  (x, y)
extent      : -3850000, 3750000, -5350000, 5850000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : nsidcCDRiceSQnhmday_935c_47bd_a147.nc 
varname     : seaice_conc_monthly_cdr (NOAA/NSIDC Climate Data Record of Passive Microwave Monthly Northern Hemisphere Sea Ice Concentration) 
name        : seaice_conc_monthly_cdr 
unit        :                       1 
time        : 2019-12-16

Как видите, номинальное разрешение ячейки 25000 x 25000.

res(r)
[1] 25000 25000

Поскольку это константа, raster не хочет вычислять ее для вас. Но terra даст вам либо наивную область (на основе постоянного номинального разрешения):

na <- cellSize(r, transform=FALSE)
minmax(na)
#         area
#[1,] 6.25e+08
#[2,] 6.25e+08

Или настоящая область. То есть с учетом искажений:

a <- cellSize(r, mask=F, unit="m", transform=TRUE)
a <- mask(a, r)
a 
#class       : SpatRaster 
#dimensions  : 448, 304, 1  (nrow, ncol, nlyr)
#resolution  : 25000, 25000  (x, y)
#extent      : -3850000, 3750000, -5350000, 5850000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :      area 
#min value   : 385538624
#max value   : 664449196 
#time        : 2019-12-16 

plot(a)

введите здесь описание изображения

Чтобы получить площадь морского льда в км2

x <- r * a / 1000000
global(x, "sum", na.rm=TRUE)
#                             sum
#_conc_monthly_cdr 11312232

Наивная оценка дает

n <- r * prod(res(r)) / 1000000
global(n, "sum", na.rm=TRUE)
#                             sum
#seaice_conc_monthly_cdr 11134025

Это примерно на 1,5% ниже

round(11134025 / 11312232, 3)
#[1] 0.984

Это не такое отличие, потому что льда нет по краям, в местах сильного искажения. А также есть некоторая компенсация в центральной области (и меньшие, и большие площади), чем 25*25 км2

person Robert Hijmans    schedule 06.05.2021
comment
Спасибо за ответ! Несмотря на то, что он говорит, что разрешение однородно, есть искажение, и ячейки неравномерны по широте. Таким образом, я считаю, что этот подход бесполезен. - person Angel Segura; 07.05.2021
comment
я расширил свой ответ - person Robert Hijmans; 07.05.2021
comment
Привет Роберт, Это потрясающе! Спасибо за расширение вашего ответа, чтобы подробно объяснить проблему. Я пройду через это обязательно. Я проголосовал за вопрос ... но пока я новичок, он не записан. С уважением - person Angel Segura; 09.05.2021