Susan O'Neil bug report

[jonathancallahan]

Hi Jon,

I have been using the raster functionality in the AirFireModeling package and noticed differences in data when I subset a RasterBrick.

Attached is ppt with a single slide with two figures. One figure is raster_leaflet of the full brick of data at a particular time, and one figure is a raster_leaflet of the subset of the brick at the same time. I’ve circled the three grid cells that do not agree. Two of those grids are fringe grids in the subset data, so I can understand them not agreeing (perhaps), but it’s the grid in the middle that made me think I should send this to you.

I’ve attached the code excerpt to recreate this (code_for_Jon.R), and put the smoke_dispersion_v2.nc at: https://haze.airfire.org/webaccess/susan/temp/forJon/

I also attached my modified raster_leaflet_smo.R so that raster_leaflet works correctly in code_for_Jon.R.

What do you think?

Thanks!

Susan

code_for_Jon.R.txt raster_leaflet_smo.R.txt ppt_for_Jon.pptx

[jonathancallahan]

The original code from Susan is for an earlier version of AirFireModeling and may need to have the following changes made:

• model ==> modelName

---

This will be an exercise in debugging.

• First, you need to see if you can reproduce the error she is seeing.
• If you can, then try to isolate precisely where the error is happening by creating the simplest reproducible example possible.
• At this point you should let Jon know what you've learned unless it is an obvious error. Those you should just fix.

[tabrasel]

This is the example Susan provided:

It does not look like there are any differences when we compare the actual cell values though:

Subset raster values (NAs present because the distance subset is circular):

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]   NA   26   33    3    0    0   NA   NA
[2,]    1   22 2292    0    0    0    0   NA
[3,]    0    0   87    5    0    0    0    0
[4,]    0    0   13    9    1    0    0   NA
[5,]   NA    0   15   93    7    5    2   NA
[6,]   NA   NA    0    0    0   NA   NA   NA

Full raster values from the same area:

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    5   26   33    3    0    0    0    0
[2,]    1   22 2292    0    0    0    0    0
[3,]    0    0   87    5    0    0    0    0
[4,]    0    0   13    9    1    0    0    0
[5,]    0    0   15   93    7    5    2    2
[6,]    0    0    0    0    0    2    0    1

The retained values match perfectly, so the coloration issue seems to stem from how the two rasters are projected and then sampled when drawn on Leaflets. If you look closely you can see that the cell grids align on the same longitude, but not on the same latitude:

The full raster:

dimensions : 134, 134, 17956  (nrow, ncol, ncell)
resolution : 0.01500001, 0.015  (x, y)
extent     : -120.0025, -117.9925, 36.9975, 39.0075  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

The subset raster:

dimensions : 6, 8, 48  (nrow, ncol, ncell)
resolution : 0.01500001, 0.015  (x, y)
extent     : -119.1775, -119.0575, 37.9125, 38.0025  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

So in order for the rasters to be perfectly aligned (and therefore be sampled and colored the same way), we might snap every raster to a predetermined grid when drawn on a Leaflet. We can use raster::shift() to translate them, but I'm not sure if this is the best solution and would work in every case.

This grid misalignment might also explain why the shape of the subset raster on a Leaflet does not exactly match the circular shape of its value matrix shown above.

[tabrasel]

Haha, this is actually way less complicated to fix. The raster_leaflet() function has been using the default "bilinear" method to color grid cells. The different raster sizes (subset or not) impact how values are interpolated between neighbors, and that seems to be the issue behind the color variations. All that needs to be done to make them match is setting the leaflet::addRasterImage() method to use "nearest neighbor" within the raster_leaflet() function. I'll make this change, but maybe we can add the method param to raster_leaflet() as well and default it to "ngb"?

This also fixes the circular shape of the raster on the Leaflet. Yay!