MazamaScience / AirFireModeling

Utilities to ease merging of USFS AirFire model output and monitoring data.
0 stars 0 forks source link

review raster grid details #39

Open jonathancallahan opened 4 years ago

jonathancallahan commented 4 years ago

Today I found out that the v2 grids we are producing are slightly off in some cases. We need to review to make sure we are getting this right.

This issue will involve communication with USFS AirFire personnel which Jon can help expedite.

The final result should be an RMarkdown report describing what you learn about model grids and their raster representation.


People involved:

Model grid information we obtained from Robert is recreated in R/AirFireModeling.R as bluesky_modelInfo. We should check with Robert that these are correct and make sure we understand how they are interpreted -- edge of cell vs center etc.

Then we should use the ncdump facility to look at grid information from the models before and after we convert to v2 to see if anything changes.


Here is the comparison of a population raster generated by Andrew Bailey using the correct grid and the grid information currently generated by AirFireModeling:

  # NOTE:  Right now we have the following mismatch.

  # NOTE:  > print(populationRaster)
  # NOTE:  class      : RasterLayer 
  # NOTE:  dimensions : 264, 514, 135696  (nrow, ncol, ncell)
  # NOTE:  resolution : 0.04, 0.04  (x, y)
  # NOTE:  extent     : -128.55, -107.99, 39.75, 50.31  (xmin, xmax, ymin, ymax)
  # NOTE:  crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
  # NOTE:  source     : /Users/jonathan/Projects/PWFSL/2020/covid-r-scripts/data/PNW-4km_pop.nc 
  # NOTE:  names      : population 
  # NOTE:  zvar       : population 

  # NOTE:  > print(modelBrick[[1]])
  # NOTE:  class      : RasterLayer 
  # NOTE:  band       : 1  (of  72  bands)
  # NOTE:  dimensions : 264, 514, 135696  (nrow, ncol, ncell)
  # NOTE:  resolution : 0.04000001, 0.03999999  (x, y)
  # NOTE:  extent     : -128.57, -108.01, 39.73, 50.29  (xmin, xmax, ymin, ymax)
  # NOTE:  crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
  # NOTE:  source     : /private/var/folders/vd/zpgw5sv92ngdx11k5dzqz5800000gn/T/RtmpJl3g6g/PNW-4km_2020061100_v2.nc 
  # NOTE:  names      : X1591837200 
  # NOTE:  z-value    : 1591837200 
  # NOTE:  zvar       : PM25 
  # NOTE:  level      : 1 

Here is a sample of how to do this work with the R "Console" and the Unix "Terminal":

> library(AirFireModeling)
> setModelDataDir('~/Data/BlueSky')
> 
> filePath <- bluesky_download(modelName = "PNW-4km", modelRun = 2019100900)
Auto-detected PNW-4km BlueSky Output Version: 1.0.0
trying URL 'https://haze.airfire.org/bluesky-daily/output/standard/PNW-4km/2019100900/forecast/data/smoke_dispersion.nc'
Content type 'application/x-netcdf' length 38545660 bytes (36.8 MB)
==================================================
downloaded 36.8 MB

> filePath
[1] "/Users/jonathan/Data/BlueSky/PNW-4km_2019100900.nc"
> v2Path <- bluesky_toCommonFormat(filePath, clean = FALSE)
> v2Path
[1] "/Users/jonathan/Data/BlueSky/PNW-4km_2019100900_v2.nc"

Now, at the Unix terminal:

$ ncdump -h /Users/jonathan/Data/BlueSky/PNW-4km_2019100900.nc
netcdf PNW-4km_2019100900 {
dimensions:
        TSTEP = UNLIMITED ; // (71 currently)
        LAY = 1 ;
        ROW = 264 ;
        COL = 514 ;
        VAR = 1 ;
        DATE-TIME = 2 ;
variables:
        float PM25(TSTEP, LAY, ROW, COL) ;
                PM25:long_name = "PM25            " ;
                PM25:units = "ug/m^3          " ;
                PM25:var_desc = "PM25                                                                            " ;
        int TFLAG(TSTEP, VAR, DATE-TIME) ;
                TFLAG:units = "<YYYYDDD,HHMMSS>" ;
                TFLAG:long_name = "TFLAG           " ;
                TFLAG:var_desc = "Timestep-valid flags:  (1) YYYYDDD or (2) HHMMSS                                " ;

// global attributes:
                :IOAPI_VERSION = "$Id: @(#) ioapi library version 3.0 $                                           " ;
                :EXEC_ID = "????????????????                                                                " ;
                :FTYPE = 1 ;
                :CDATE = 2019282 ;
                :CTIME = 92601 ;
                :WDATE = 2019282 ;
                :WTIME = 92601 ;
                :SDATE = 2019282 ;
                :STIME = 10000 ;
                :TSTEP = 10000 ;
                :NTHIK = 1 ;
                :NCOLS = 514 ;
                :NROWS = 264 ;
                :NLAYS = 1 ;
                :NVARS = 1 ;
                :GDTYP = 1 ;
                :P_ALP = 0. ;
                :P_BET = 0. ;
                :P_GAM = 0. ;
                :XCENT = -118.290000915527 ;
                :YCENT = 45.0099983215332 ;
                :XORIG = -128.550003051758 ;
                :YORIG = 39.75 ;
                :XCELL = 0.0399999991059303 ;
                :YCELL = 0.0399999991059303 ;
                :VGTYP = 5 ;
                :VGTOP = -9999.f ;
                :VGLVLS = 100.f, 0.f ;
                :GDNAM = "HYSPLIT CONC    " ;
                :UPNAM = "hysplit2netCDF  " ;
                :VAR-LIST = "PM25            " ;
                :FILEDESC = "Hysplit Concentration Model Output                                              lat-lon coordinate system
...

and, to see all the coordinate values:

ncdump -c /Users/jonathan/Data/BlueSky/PNW-4km_2019100900_v2.nc
...

So far it all looks good.

Reading things in with raster::rasterBrick() gives us the suspect grid info:

> rasterBrick <- raster::brick(v2Path)
> print(rasterBrick)
class      : RasterBrick 
dimensions : 264, 514, 135696, 71  (nrow, ncol, ncell, nlayers)
resolution : 0.04000001, 0.03999999  (x, y)
extent     : -128.57, -108.01, 39.73, 50.29  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : /Users/jonathan/Data/BlueSky/PNW-4km_2019100900_v2.nc 
names      : X1570582800, X1570586400, X1570590000, X1570593600, X1570597200, X1570600800, X1570604400, X1570608000, X1570611600, X1570615200, X1570618800, X1570622400, X1570626000, X1570629600, X1570633200, ... 
z-value    : 1570582800, 1570834800 (min, max)
varname    : PM25 
level      : 1 
...
> raster::extent(rasterBrick)
class      : Extent 
xmin       : -128.57 
xmax       : -108.01 
ymin       : 39.73 
ymax       : 50.29 

Why is xmin all of a sudden -128.57 instead of -128.55?

tabrasel commented 4 years ago

From what I've picked up so far, I think that our grids are actually defined correctly. I'm pretty sure that--like you said--the issue might stem from cells being placed at their centers rather than corners. Unfortunately, the link to Andrew's grid is down right now so I can't tell yet why there's a difference between his and ours. I'll try to explain my thinking though using the NetCDF you demoed (modelName = "PNW-4km", modelRun = 2019100900):

We can find the NetCDF's origin (south-westernmost coordinate) in the original NetCDF with ncdump -h ~/Data/BlueSky/PNW-4km_2019100900.nc:

## netcdf PNW-4km_2019100900 {
##      ...
##      :XORIG = -128.550003051758 ;
##      :YORIG = 39.75 ;
##.     ...

Let's convert the NetCDF, load it as a raster, and check out its extent:

> filePath <- bluesky_download(modelName = "PNW-4km", modelRun = 2019100900)
> v2Path <- bluesky_toCommonFormat(filePath, clean = FALSE)
> raster <- bluesky_load(localPath = v2Path)
> raster::extent(raster)
class      : Extent 
xmin       : -128.57 
xmax       : -108.01 
ymin       : 39.73 
ymax       : 50.289997 

The extent's south-westernmost coordinate (-128.57, 39.73) is different from the original NetCDF's (-128.550003051758, 39.75). It's shifted by (-0.01999694824, -0.02) which, incidentally, is about half of the raster's cell resolution (0.040000008, 0.039999987). That doesn't mean that anything is wrong with our grid placement though. Let's look at the actual coordinates of the South-westernmost cell:

> # South-westernmost cell's long/lat coordinates
> sw <- raster::cellFromRowCol(raster, raster@nrows, 1)
> raster::xyFromCell(raster, sw)
           x         y
[1,] -128.55 39.75

The same as the original NetCDF! So conversion doesn't seem to be mutating the actual coordinates. So why was the extent of the raster off? Well, let's look at the extent of just the south-westermost cell:

> raster::extentFromCells(raster, sw)
class      : Extent 
xmin       : -128.57 
xmax       : -128.53 
ymin       : 39.73 
ymax       : 39.77 

The extent of a single cell is offset by the half the cell's resolution on each side. Therefore, cells are centered on their coordinates, and that's what the case is with the extent of the whole raster. It's borders surround the cell coordinates it contains with a bit of padding:

+---------------+
| *   *   *   * |
|               |
| *   *   *   * |
|               |
| *   *   *   * |
+---------------+

So unless our actual NetCDF coordinates are different from Andrew's, it seems that our raster is drawing its cells where it should be. I'm a bit confused then how Andrew's extent could contain (-128.55, 39.75) when cells are supposed to be defined at their center, not at their corner.