使用立体极坐标投影从卫星数据估计网格单元面积 一些 R 试验

编程问答 2022-05-03

问题描述

我正在尝试处理来自极地地区的卫星数据。它们可以下载为 .nc（netcdf，见下文）。 它们是在规则网格 (https://nsidc.org/data/polar-stereo/ps_grids.html) 中的极坐标立体投影。 我想通过将每个单元格中冰盖的分数乘以单元格的面积来估计每个单元格的面积以计算冰盖的面积。

只有当统一网格处于地理坐标（经纬度）时，我才能使用 {raster} 中的“area”提取单元格区域。但是，当网格在立体坐标中是均匀的时，我找不到任何 R 函数来执行此操作。

我对空间数据和投影非常陌生，也许我遗漏了一些重要的东西。
下面是一些相关信息和一段带有两次试验的代码。

相关数据和信息来源

极地手表 https://polarwatch.noaa.gov/erddap/

海冰浓度数据 下载数据的网址

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)]

可以使用 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 公里 数据集 ID：nsidcCDRice_sh_grid 网址 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)]

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

一些 R 试验

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

任何帮助将不胜感激。

提前致谢 天使

解决方法

这里我用 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)
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 <- area(r,correct=FALSE,sum=FALSE)
minmax(na)
#         area
#[1,] 6.25e+08
#[2,] 6.25e+08

或者真实区域。也就是说，考虑失真：

a <- area(r,correct=TRUE,mask=TRUE,sum=FALSE)
a 
#class       : SpatRaster 
#dimensions  : 448,nlyr)
#resolution  : 25000,y)
#extent      : -3850000,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 

plot(a)

enter image description here

以km为单位获取海冰面积2

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

朴素的估计收益

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

大约低了 1.5%

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

这与那没有什么不同，因为在严重变形的区域边缘没有冰。并且在中心区（小区和大区）也有一些补偿，超过25*25km2

areaareapolar-coordinatesrrsatellite