使用立体极坐标投影从卫星数据估计网格单元面积 [英] Estimate grid cell area from satelite data using a Stereographic polar projection

问题描述

我正在尝试处理来自极地地区的卫星数据.它们可以下载为 .nc(netcdf，见下文).它们在规则网格中的极地立体投影(

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

x <- r * a/1000000全局(x，总和"，na.rm=TRUE)#总和#_conc_monthly_cdr 11312232

朴素的估计收益

n <- r * prod(res(r))/1000000全局(n，总和"，na.rm=TRUE)#总和#seaice_conc_monthly_cdr 11134025

大约低了 1.5%

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

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

I am trying to work with satelite data from polar regions. They can be download as .nc (netcdf, see below). They are in polar stereographic projection in a regular grid (https://nsidc.org/data/polar-stereo/ps_grids.html). I would like to estimate each cell area to calculate the area of ice cover by multiplying the fraction of ice cover in each cell with the area of the cell.

I am able to extract cell area using "area" from {raster} only when the uniform grid is in geographical coordinates (lat-lon). However I could not find any R function to do this when the grid is uniform in stereographic coordinates.

I am very new to spatial data and projections and perhaps I am missing something important.
Below some relevant information and a piece of code with two trials.

RELEVANT DATA AND INFORMATION SOURCES

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

Sea Ice concentration data url to downoald data

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

Can be download in R using

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

Data for the grid Sea Ice Concentration Lat-Lon Grid, NOAA/NSIDC Climate Data Record V3, Antarctic, 25km Dataset ID: nsidcCDRice_sh_grid url 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 Table 4. Southern Hemisphere Projection Based on WGS 1984 Latitude of True Origin -70

SOME R TRIALS

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)

which gives

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

Using the grid provided in the data sources

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

But gives

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

Any help will be appreciated.

Thanks in advance Angel

解决方案

Here I use terra the replacement for 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 

As you can see, the nominal cell resolution is 25000 x 25000

res(r)
[1] 25000 25000

Because it is constant, raster does not want to compute it for you. But terra will give you either the naïve area (based on the constant the nominal resolution):

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

Or the the true area. That is, accounting for distortion:

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)

To get the sea-ice area in km²

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

The naïve estimate yields

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

That is about 1.5% lower

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

This is not that different, because there is no ice at the edges, in the areas of severe distortion. And there is also some compensation in the central area (both smaller and larger areas) than 25*25 km²

这篇关于使用立体极坐标投影从卫星数据估计网格单元面积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！