Problem with "striping" when putting rectangular grids into geom_sf()

[rmendels]

The following problem was noticed when working the "plotdap" package that is under development at:

https://github.com/ropensci/plotdap

To recreate the problem, you will need the package rerddap as well as the plotdap package. I have also attached the output from using geom_sf() and using base graphics, two options in plotdap.

Try this code:

library(rerddap)
library(plotdap)
library(ggplot2)
sstInfo <- info('jplMURSST41')
murSST <- griddap(sstInfo, latitude = c(22., 51.), longitude = c(-140., -105), time = c('last','last'), fields = 'analysed_sst')
add_griddap(plotdap(), murSST, ~analysed_sst, maxpixels = 50000)

Note that there is a pronounced "striping" in both directions. You can generate the plot using base graphics by changing the last lie to:

add_griddap(plotdap("base"), murSST, ~analysed_sst, maxpixels = 50000)

Compare also with (though the color gradient is different):

ggplot(data = murSST$data, aes(x = lon, y = lat, fill = analysed_sst)) + 
    geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
    geom_raster(interpolate = FALSE, na.rm = TRUE) + theme_bw() + ylab("latitude") + xlab("longitude") +
     coord_fixed(1.3, xlim = c(-140, -105),  ylim = c(22., 51.)) + ggtitle("Latest MUR SST")

Notice the difference there also. I brought this issue up with Carson Sievert who developed plotdap, and he commented that the problem was that it "It's an unfortunate artifact of using geom_sf() to draw "square" rasters via filled polygons"

The combination of rerddap and pltodap makes it very easy for working scientists who don't program much to quickly download and map data, and I prefer ggplot2 for my graphics, for the flexibility it allows after the initial generation (plotdap allows you to modify the graphic as usual in ggplot2). It would be great if a solution can be found for this.

Thanks,

-Roy

[edzer]

Thanks! Strange enough, I see

with:

``` > devtools::session_info() Session info ------------------------------------------------------------------- setting value version R version 3.4.1 (2017-06-30) system x86_64, linux-gnu ui X11 language en_US collate en_US.UTF-8 tz Europe/Berlin date 2017-07-24 Packages ----------------------------------------------------------------------- package * version date source assertthat 0.2.0 2017-04-11 CRAN (R 3.4.0) bindr 0.1 2016-11-13 cran (@0.1) bindrcpp 0.2 2017-06-17 CRAN (R 3.4.0) colorspace 1.3-2 2016-12-14 CRAN (R 3.4.0) curl 2.8.1 2017-07-21 cran (@2.8.1) data.table 1.10.4 2017-02-01 CRAN (R 3.4.0) DBI 0.7 2017-06-18 CRAN (R 3.4.0) devtools 1.12.0 2016-12-05 CRAN (R 3.4.0) digest 0.6.12 2017-01-27 CRAN (R 3.4.0) dplyr 0.7.2 2017-07-20 cran (@0.7.2) foreign 0.8-69 2017-06-21 CRAN (R 3.4.0) ggplot2 * 2.2.1.9000 2017-07-24 Github (tidyverse/ggplot2@331977e) glue 1.1.1 2017-06-21 CRAN (R 3.4.0) gtable 0.2.0 2016-02-26 CRAN (R 3.4.0) hoardr 0.2.0 2017-05-10 cran (@0.2.0) httr 1.2.1 2016-07-03 CRAN (R 3.4.0) jsonlite 1.5 2017-06-01 cran (@1.5) labeling 0.3 2014-08-23 CRAN (R 3.4.0) lattice 0.20-35 2017-03-25 CRAN (R 3.3.3) lazyeval 0.2.0 2016-06-12 CRAN (R 3.4.0) magrittr 1.5 2014-11-22 CRAN (R 3.4.0) mapdata * 2.2-6 2016-01-14 CRAN (R 3.4.0) maps * 3.2.0 2017-06-08 cran (@3.2.0) maptools 0.9-2 2017-03-25 CRAN (R 3.4.0) memoise 1.1.0 2017-04-21 CRAN (R 3.4.0) munsell 0.4.3 2016-02-13 CRAN (R 3.4.0) ncdf4 1.16 2017-04-01 CRAN (R 3.4.0) pkgconfig 2.0.1 2017-03-21 cran (@2.0.1) plotdap * 0.0.1 2017-07-24 Github (ropensci/plotdap@3f167c1) plyr 1.8.4 2016-06-08 CRAN (R 3.4.0) R6 2.2.2 2017-06-17 cran (@2.2.2) rappdirs 0.3.1 2016-03-28 CRAN (R 3.4.0) raster 2.5-8 2016-06-02 CRAN (R 3.4.0) Rcpp 0.12.12 2017-07-15 CRAN (R 3.4.1) rerddap * 0.4.2 2017-05-12 cran (@0.4.2) rgdal 1.2-8 2017-07-01 CRAN (R 3.4.1) rgeos 0.3-23 2017-04-06 CRAN (R 3.4.0) rlang 0.1.1.9000 2017-07-24 Github (hadley/rlang@342c473) scales 0.4.1.9002 2017-07-24 Github (hadley/scales@6db7b6f) sf * 0.5-2 2017-07-12 CRAN (R 3.4.1) sp 1.2-5 2017-06-29 CRAN (R 3.4.1) tibble 1.3.3 2017-05-28 CRAN (R 3.4.0) tidyr 0.6.3 2017-05-15 cran (@0.6.3) udunits2 0.13 2016-11-17 CRAN (R 3.4.0) units 0.4-6 2017-07-12 local withr 1.0.2 2016-06-20 CRAN (R 3.4.0) xml2 1.1.1 2017-01-24 CRAN (R 3.4.0) ```

With sf from github it breaks; will look into why.

[edzer]

@mdsumner is isolating the raster gg plotting still on your todo list (mentioned here), or has it been resolved?

[rmendels]

Thanks. It has been awhile since I first noticed this (mid-June). I have been meaning to raise this issue, but just got around to it, maybe it has been resolved. I only recently figured out how to install SF from source, so that I could use the development version, and don't always update ggplot2 from development. Can I ask what version of everything you are using (which version of SF and which version of ggplot2)? I will try to install them and do more testing.

Thanks for looking into this. The image you sent is much more like we are after.

-Roy

[edzer]

For versions, click on the issue, and then on the "details" section in my first answer.

[edzer]

The problem I saw was in my fork of ggplot2; you won't need that.

[cpsievert]

I get the same issue as @rmendels with a similar setup:

```r Session info ------------------------------------------------------------------ setting value version R version 3.4.1 (2017-06-30) system x86_64, darwin15.6.0 ui X11 language (EN) collate en_US.UTF-8 tz America/Chicago date 2017-07-24 Packages ---------------------------------------------------------------------- package * version date source assertthat 0.2.0 2017-04-11 CRAN (R 3.4.0) autoinst 0.0.0.9000 2017-04-27 Github (jimhester/autoinst@dfbdc50) base * 3.4.1 2017-07-07 local bindr 0.1 2016-11-13 CRAN (R 3.4.0) bindrcpp 0.2 2017-06-17 CRAN (R 3.4.0) colorout 1.1-2 2017-04-23 Github (jalvesaq/colorout@020a14d) colorspace 1.3-2 2016-12-14 CRAN (R 3.4.0) compiler 3.4.1 2017-07-07 local curl 2.8.1 2017-07-21 CRAN (R 3.4.1) data.table 1.10.4 2017-02-01 CRAN (R 3.4.0) datasets * 3.4.1 2017-07-07 local DBI 0.7 2017-06-18 CRAN (R 3.4.0) devtools * 1.13.2 2017-06-02 CRAN (R 3.4.0) digest 0.6.12 2017-01-27 CRAN (R 3.4.0) dplyr 0.7.2 2017-07-20 CRAN (R 3.4.1) foreign 0.8-69 2017-06-22 CRAN (R 3.4.1) fortunes 1.5-4 2016-12-29 CRAN (R 3.4.0) ggplot2 * 2.2.1.9000 2017-07-24 Github (hadley/ggplot2@331977e) git2r 0.19.0 2017-07-19 CRAN (R 3.4.1) glue 1.1.1 2017-06-21 CRAN (R 3.4.0) graphics * 3.4.1 2017-07-07 local grDevices * 3.4.1 2017-07-07 local grid 3.4.1 2017-07-07 local gtable 0.2.0 2016-02-26 CRAN (R 3.4.0) hoardr 0.2.0 2017-05-10 CRAN (R 3.4.0) htmltools 0.3.6 2017-04-28 CRAN (R 3.4.0) htmlwidgets 0.9 2017-06-06 Github (ramnathv/htmlwidgets@387508c) httr 1.2.1 2016-07-03 CRAN (R 3.4.0) inline 0.3.14 2015-04-13 CRAN (R 3.4.0) jsonlite 1.5 2017-06-01 CRAN (R 3.4.0) knitr 1.16 2017-05-18 CRAN (R 3.4.0) labeling 0.3 2014-08-23 CRAN (R 3.4.0) lattice 0.20-35 2017-03-25 CRAN (R 3.4.1) lazyeval 0.2.0.9000 2017-05-08 Github (hadley/lazyeval@c155c3d) magrittr 1.5 2014-11-22 CRAN (R 3.4.0) mapdata * 2.2-6 2016-01-14 CRAN (R 3.4.0) maps * 3.2.0 2017-06-08 CRAN (R 3.4.0) maptools 0.9-2 2017-03-25 CRAN (R 3.4.0) memoise 1.1.0 2017-04-21 CRAN (R 3.4.0) methods * 3.4.1 2017-07-07 local munsell 0.4.3 2016-02-13 CRAN (R 3.4.0) ncdf4 1.16 2017-04-01 CRAN (R 3.4.0) pkgconfig 2.0.1 2017-03-21 CRAN (R 3.4.0) pkgload 0.0.0.9000 2017-04-23 Github (r-pkgs/pkgload@9093b96) plotdap * 0.0.1 2017-07-24 Github (ropensci/plotdap@3f167c1) plyr 1.8.4 2016-06-08 CRAN (R 3.4.0) R6 2.2.2 2017-06-17 CRAN (R 3.4.0) rappdirs 0.3.1 2016-03-28 CRAN (R 3.4.0) raster 2.5-8 2016-06-02 CRAN (R 3.4.0) Rcpp 0.12.12 2017-07-15 CRAN (R 3.4.1) reprex * 0.1.1 2017-01-13 CRAN (R 3.4.0) rerddap * 0.4.2.9130 2017-07-24 Github (ropensci/rerddap@6d7ddad) rgdal 1.2-8 2017-07-01 CRAN (R 3.4.1) rgeos 0.3-23 2017-04-06 CRAN (R 3.4.0) rlang 0.1.1 2017-05-18 CRAN (R 3.4.0) scales 0.4.1.9002 2017-07-17 Github (hadley/scales@6db7b6f) sf 0.5-2 2017-07-12 CRAN (R 3.4.1) sp 1.2-5 2017-06-29 CRAN (R 3.4.1) stats * 3.4.1 2017-07-07 local tibble 1.3.3 2017-05-28 CRAN (R 3.4.0) tidyr 0.6.3 2017-05-15 CRAN (R 3.4.0) tools 3.4.1 2017-07-07 local udunits2 0.13 2016-11-17 CRAN (R 3.4.0) units 0.4-5 2017-06-15 CRAN (R 3.4.0) utils * 3.4.1 2017-07-07 local withr 1.0.2 2016-06-20 CRAN (R 3.4.0) xml2 1.1.1 2017-01-24 CRAN (R 3.4.0) ```

[rmendels]

I run on a Mac, and I believe Carson does also. I believe Edzer runs on Linux, is that correct? Could this be a OS related issue?

Just a thought. I haven't had a chance to check that I have all the right versions and rerun.

-Roy

[rmendels]

I just ran on my Mac, I have the same versions of ggplot2 and sf. I still get the striping. I thought it might be related to a problem that geom_raster() had if the lat-lon wasn't precisely even, but I tried a dataset I know has even lat-lon grids, and that shows it too. This one is not as fine-grained, it what it looks like is the outlines of each cell are being drawn in a given color, or maybe it i each lat-lon being drawn. Try this:

library(rerddap) library(plotdap) library(ggplot2) library(sf) dataInfo <- rerddap::info('hawaii_d90f_20ee_c4cb') xpos <- c(210.25, 240.25) ypos <- c(20.25, 60.25) zpos <- c(70.02, 70.02) tpos <- c('2010-12-15', '2010-12-15') soda70 <- griddap(dataInfo, longitude = xpos, latitude = ypos, time = tpos, depth = zpos, fields = 'temp' ) plotdap() %>% add_griddap(soda70, ~temp, maxpixels = 50000)

The graphic I get is attached, if a Linux box get different that will tell us something. This dataset is on a 0.5 degree grid, the other were on like a 0.01 degree grid.

Thanks,

-Roy

[mdsumner]

@edzer I don't want to prolong this issue, I don't see it's relevant to sf - but here's my thoughts. I'll be exploring this example a lot more. Getting an efficient model for general gridded data for the tidyverse is very much an active concern, I'm exploring this generally here: https://github.com/hypertidy - gridcol is the most recent attempt at an idea for a table-raster, but it's only suitable for affine rasters - ultimately a multi-table approach is needed to deal with relational data, which geotransforms and rectlinear axes and curvlinear grids are forms of.

For this graphic in particular, it might be relevant that the coordinate arrays are truly "rectlinear" in a very minute sense, though not enough for us to really notice.

range(diff(sort(unique(murSST$data$lon))))
[1] 0.009994507 0.010009766
> range(diff(sort(unique(murSST$data$lat))))
[1] 0.009998322 0.010002136

I am interested to explore this generally, it's very much a dark corner in R that's not widely discussed. Geom_raster and friends is somehow hiding this non-regular-affine property, but finding out where it gets handled or ignored would be interesting for my purposes. Perhaps it's just passing down to geom_tile, perhaps it swallows the tiny numeric differences silently just as raster does.

Certainly raster ignores the issue totally, and doesn't detect or care about the fine discrepancy from regularity here:

f <- attr(murSST, "path")
library(raster)
## note no complaint about "non-regular axes" here
d <- raster(f)
d
class       : RasterLayer 
dimensions  : 2901, 3501, 10156401  (nrow, ncol, ncell)
resolution  : 0.01, 0.01  (x, y)
extent      : -140.005, -104.995, 21.995, 51.005  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

## confirm we see the same non-regular property direct from the source
conn <- ncdf4::nc_open(f)
conn <- ncdf4::nc_open(f)
range(diff(ncdf4::ncvar_get(conn, "longitude")))
[1] 0.009994507 0.010009766
range(diff(ncdf4::ncvar_get(conn, "latitude")))
[1] 0.009998322 0.010002136

I also don't see how this is related to geom_sf - but happy to be corrected, I'd expect the same striping to be a product of geom_raster and a specific graphics device with a simpler reprex from an anologous data set. I assume the problem is seen using the quartz device? What about png(), pdf() on MacOS? I'd first see if I can reproduce from a simple data set with a similar fine "irregularity" in the axes, and whether it is resolved by "orthogonalizing" the coordinate values to the intended 0.01 grain.

[rmendels]

Hi Michael;

Thank you for the comments. I thought about the slight non-regularity also, as I had run into that previously with geom_raster(). Please look at the soda70 example I sent later. That should be on an absolutely regular grid, and I still get the striping. Even more, and this is way outside of anything I know about, is that I don't see the same thing with base graphics. I won't get to this tonight, but I will take one of the examples, and redo the lats and lons stored in the file to be generated so therefore hopefully totally regular, and see whether that makes the lines go away. The other thing I am wondering is rerddap actually gets from the service a netcdf file, and then reads it in and "melts" it to a long form. I have wondered whether that affects things

The other thing that is confusing is the third part I sent, where I left out a line that is needed. It is the example that uses geom_raster() directly and does not produce the striping. The complete code would be:

library(rerddap) library(ggplot2) library(mapdata) w <- map_data("worldHires", ylim = c(22., 51.), xlim = c(-140, -105)) sstInfo <- info('jplMURSST41') murSST <- griddap(sstInfo, latitude = c(22., 51.), longitude = c(-140., -105), time = c('last','last'), fields = 'analysed_sst') ggplot(data = murSST$data, aes(x = lon, y = lat, fill = analysed_sst)) + geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + geom_raster(interpolate = FALSE, na.rm = TRUE) + theme_bw() + ylab("latitude") + xlab("longitude") + coord_fixed(1.3, xlim = c(-140, -105), ylim = c(22., 51.)) + ggtitle("Latest MUR SST")

The color bar is different, but you will notice there is no striping (picture below)

Whether this is an SF problem or a ggplot2 issue, I also do not know. I started with Edzer as I had communicated with him before, and also I wasn't sure who was doing the development on geom_sf. The type of services and data we are getting are very common in meteorology and oceanography. ERDDAP servers alone allow you to access well over a petabyte of data. We have shown plotdap in demos to a number of people who were very excited about it. We would like it to produce as high as quality graphics as possible.

So again, I thank everyone for their help on this. The development effort is much appreciated.

-Roy

[edzer]

I believe that the original problem here involves converting grid cells to POLYGON, then plotting them through sf (using, in the end grid::pathGrob); the regular gridded one I assume uses in the end grid::rasterGrob.

@rmendels the graphics you mention didn't end up on the GH issue: did you maybe attach them to an email?

[mdsumner]

Ah, I see. Sorry for the noise, but thanks for the example, it's interesting for my purposes I'll summarize and report independently. (FWIW, raster does complain about the irregular spacing, it's just that the download path file above has no coordinates with it, or are missed by raster. We have the murSST on hand if anyone wants more direct access to it via RStudio or ssh, the total size - since 2002, daily - is 2000Gb. A simple script can get it all or any part if you wanted it locally: https://github.com/AustralianAntarcticDivision/bowerbird/blob/master/README.md#ghrsst-level-4-mur-global-foundation-sst-v41)

[rmendels]

yes I attached them on an email. I will get into the issues and post them on the web page. Sorry, I thought if they were in an email they would make it to the web page. I need to catch up tho the thread since I went to bed last night.

Thanks again to everyone for the insight.

-Roy

On Jul 25, 2017, at 12:01 AM, Edzer Pebesma notifications@github.com wrote:

I believe that the original problem here involves converting grid cells to POLYGON, then plotting them through sf (using, in the end grid::pathGrob); the regular gridded one I assume uses in the end grid::rasterGrob.

@rmendels the graphics you mention didn't end up on the GH issue: did you maybe attach them to an email?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

---

"The contents of this message do not reflect any position of the U.S. Government or NOAA."

---

Roy Mendelssohn Supervisory Operations Research Analyst NOAA/NMFS Environmental Research Division Southwest Fisheries Science Center Note new street address 110 McAllister Way Santa Cruz, CA 95060 Phone: (831)-420-3666 Fax: (831) 420-3980 e-mail: Roy.Mendelssohn@noaa.gov www: http://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill." "From those who have been given much, much will be expected" "the arc of the moral universe is long, but it bends toward justice" -MLK Jr.

[rmendels]

The images I had attached on the emails didn't make it.

The following code is on a gird that is not irregular at all:

library(rerddap) library(plotdap) library(ggplot2) library(sf) dataInfo <- rerddap::info('hawaii_d90f_20ee_c4cb') xpos <- c(210.25, 240.25) ypos <- c(20.25, 60.25) zpos <- c(70.02, 70.02) tpos <- c('2010-12-15', '2010-12-15') soda70 <- griddap(dataInfo, longitude = xpos, latitude = ypos, time = tpos, depth = zpos, fields = 'temp' ) plotdap() %>% add_griddap(soda70, ~temp, maxpixels = 50000)

It produces the following plot:

Michael Sumner asked what happens if I use the following:

plotdap() %>% add_griddap(soda70, ~temp, maxpixels = 50000, colour = "transparent")

I get:

I also pointed out that even for the MUR data if I use geom_raster() directly I do not get the striping (here the color palette will be different):

library(rerddap) library(ggplot2) library(mapdata) w <- map_data("worldHires", ylim = c(22., 51.), xlim = c(-140, -105)) sstInfo <- info('jplMURSST41') murSST <- griddap(sstInfo, latitude = c(22., 51.), longitude = c(-140., -105), time = c('last','last'), fields = 'analysed_sst') ggplot(data = murSST$data, aes(x = lon, y = lat, fill = analysed_sst)) + geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + geom_raster(interpolate = FALSE, na.rm = TRUE) + theme_bw() + ylab("latitude") + xlab("longitude") + coord_fixed(1.3, xlim = c(-140, -105), ylim = c(22., 51.)) + ggtitle("Latest MUR SST")

[cpsievert]

I believe that the original problem here involves converting grid cells to POLYGON, then plotting them through sf (using, in the end grid::pathGrob); the regular gridded one I assume uses in the end grid::rasterGrob.

Yes @edzer, I had/have the same hunch. Here is a more minimal example of the same problem. Do you know of a better approach for plotting rasters alongside other sf feature layers (via geom_sf())?

library(raster)
library(sf)
library(ggplot2)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- st_as_sf(rasterToPolygons(r))
ggplot() + geom_sf(data = s, aes(fill = test), colour = "transparent")

[rmendels]

On Jul 25, 2017, at 12:52 AM, Michael Sumner notifications@github.com wrote:

We have the murSST on hand if anyone wants more direct access to it via RStudio or ssh, the total size - since 2002, daily - is 2000Gb. A simple script can get it all or any part if you wanted it locally:

rerddap allows you to do the same. Moreover, it has a cache. and if you look in the cache there is a netcdf file stored, so you can just copy it somewhere. Moreover, if you make the same request, it doesn't actually get the data again.

-Roy

---

"The contents of this message do not reflect any position of the U.S. Government or NOAA."

---

Roy Mendelssohn Supervisory Operations Research Analyst NOAA/NMFS Environmental Research Division Southwest Fisheries Science Center Note new street address 110 McAllister Way Santa Cruz, CA 95060 Phone: (831)-420-3666 Fax: (831) 420-3980 e-mail: Roy.Mendelssohn@noaa.gov www: http://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill." "From those who have been given much, much will be expected" "the arc of the moral universe is long, but it bends toward justice" -MLK Jr.

[tim-salabim]

I thought it had to do with the size of the plotting device (as I had issues with this before), but I don't seem to have any issues with @cpsievert 's example (at any device size):

R version 3.4.1 (2017-06-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.2 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
 [1] LC_CTYPE=en_US.UTF-8      
 [2] LC_NUMERIC=C              
 [3] LC_TIME=de_DE.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=de_DE.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=de_DE.UTF-8      
 [8] LC_NAME=C                 
 [9] LC_ADDRESS=C              
[10] LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_DE.UTF-8
[12] LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets 
[6] methods   base     

other attached packages:
[1] ggplot2_2.2.1.9000 sf_0.5-2          
[3] raster_2.5-8       sp_1.2-5          

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12      lattice_0.20-35  
 [3] digest_0.6.12     plyr_1.8.4       
 [5] grid_3.4.1        gtable_0.2.0     
 [7] DBI_0.7           magrittr_1.5     
 [9] units_0.4-5       scales_0.4.1.9002
[11] rlang_0.1.1       lazyeval_0.2.0   
[13] labeling_0.3      rgdal_1.2-7      
[15] tools_3.4.1       udunits2_0.13    
[17] munsell_0.4.3     yaml_2.1.14      
[19] compiler_3.4.1    colorspace_1.3-2 
[21] tibble_1.3.3    

[rmendels]

Oh my, this appears to be a consistent pattern - no problem on Linux, a problem on the Mac. Does anyone have a Windows machine that they can test on? It would be nice to isolate where and when this happens.

-Roy

[mdsumner]

geom_raster is what I would use:

cl <- sf::st_as_sf(rasterToContour(r, level = 1000))
#s <- st_as_sf(rasterToPolygons(r))
ggplot() + 
  geom_raster(data = as.data.frame(r, xy = TRUE), aes(x, y, fill = test)) + 
  geom_sf(data = cl, colour = "white") 

I'm still exploring plotdap for the overall approach, which is very good learning for me!

My meta-answer also is "don't convert a raster to polygons (unless it's really necessary)". If you really need to there are more efficient ways, one is spex::polygonize - but still I think it's the wrong approach since you have 5 coordinates for each pixel in sf, and only 1 if geom_raster pushes direct to rasterGrob as designed.. If a different projection is required, I see that plotdap is already transforming the raster anyway so that's not the reason for conversion to polygons - is that right? Under realistic transformation the polygons will require densification e.g. with st_segmentize but even then it won't align to loxodromes without extra care and gaps can creep in as the boundaries are not shared.

One way to reduce the detail in this (immensely) high-res data is with raster's cell abstraction logic and standard data frame summarize. Data frames as handlers for rasters can be very powerful, but it needs special care and it's not at all well supported in R generally.

maxpixels <- 50000
dres <- c(mean(diff(sort(unique(murSST$data$lon)))), mean(diff(sort(unique(murSST$data$lat)))))
library(raster)
r <- raster(extent(range(murSST$data$lon) + c(-1, 1) * dres[1]/2, range(murSST$data$lat) + c(-1, 1) * dres[2]/2),
       res = dres, crs = "+init=epsg:4326")

dim(r) <- dim(r)[1:2] %/% sqrt(ceiling(ncell(r) / maxpixels))

dat <- murSST$data %>%
  mutate(bigcell = cellFromXY(r, cbind(lon, lat))) %>%
  group_by(time, bigcell) %>%
  summarize(analysed_sst = mean(analysed_sst, na.rm = FALSE)) %>%
  ungroup() %>%
  mutate(lon = xFromCell(r, bigcell), lat = yFromCell(r, bigcell))

(See how we need to juggle the r object as specification of grid topology, it's not at all easy to do this in a single-data frame context, probably not possible efficiently).

(thanks @rmendels I see all those features, it's very good - I get a bit sidetracked on these topics as I end up making connections in lots of different directions, I am writing this one up in detail now - I was definitely off-track with the task at hand)

[cpsievert]

If a different projection is required, I see that plotdap is already transforming the raster anyway so that's not the reason for conversion to polygons - is that right?

Yes, it does -- the reason for converting raster to sf is because geom_raster() doesn't support non-cartesian coordinate systems

[rmendels]

One last comment. I updated my Ubuntu that I run as a virtual machine on my Mac, and got everything installed (with some effort). I could run the commands side by side. Every example above had the same result - on the Mac, striping, on the Ubuntu running as a VM on the mac, no striping.

This is way above anything I know about.

[mdsumner]

@cpsievert I understand the reprojection problem with geom_raster, I'm saying the way to do that best is to use the facilities designed for raster - but also apply some tricks to make the transition between array and data frame forms a little easier for developers. This is the space I generally want to see more exploration in, and as ever I'm happy to do that with anyone.

https://gist.github.com/mdsumner/573ec70014df177baa2d1df7e55e1943

plotdap is already using projectRaster so the extra conversion to polygons is not required, but control over the detail-reduction certainly is. If you are interested I'll craft a PR for plotdap along those lines. The resampling applied by projectRaster is by bilinear interpolation, and it may be better to control that more explicitly to allow aggregate and other methods instead, similar to the piped example above but with an extra inverse coord transformation.

@rmendels I strongly sympathize, it's an absolute maze this stuff - at the moment it's really best to learn the available componentry rather than rely on these higher level tools. There's no language for getting between ggplot2 and spatial forms, but there are very efficient and high fidelity ways to go about it.

[rmendels]

But why does everything work okay on Ubuntu? I still wish someone could test all this on Windows. I would be very curious if this is a Mac specific issue.

Otherwise, thanks for the help.

-Roy

On Jul 25, 2017, at 5:03 PM, Michael Sumner notifications@github.com wrote:

@cpsievert I understand the reprojection problem with geom_raster, I'm saying the way to do that best is to use the facilities designed for raster - but also apply some tricks to make the transition between array and data frame forms a little easier for developers. This is the space I generally want to see more exploration in, and as ever I'm happy to do that with anyone.

https://gist.github.com/mdsumner/573ec70014df177baa2d1df7e55e1943

plotdap is already using projectRaster so I the extra conversion to polygons is not required. If you are interested I'll craft a PR for plotdap along those lines. The resampling applied by projectRaster is by bilinear interpolation, and it may be better to control that more explicitly to allow aggregate and other methods instead, similar to the piped example above but with an extra inverse coord transformation.

@rmendels I strongly sympathize, it's an absolute maze this stuff - at the moment it's really best to learn the available componentry rather than rely on these higher level tools. There's no language for getting between ggplot2 and spatial forms, but there are very efficient and high fidelity ways to go about it.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

---

"The contents of this message do not reflect any position of the U.S. Government or NOAA."

---

Roy Mendelssohn Supervisory Operations Research Analyst NOAA/NMFS Environmental Research Division Southwest Fisheries Science Center Note new street address 110 McAllister Way Santa Cruz, CA 95060 Phone: (831)-420-3666 Fax: (831) 420-3980 e-mail: Roy.Mendelssohn@noaa.gov www: http://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill." "From those who have been given much, much will be expected" "the arc of the moral universe is long, but it bends toward justice" -MLK Jr.

[rmendels]

@mdsumner Just looked at the Gist. No striping on the mac. But what is the gray area to the left of the image? That should not be there.

On Jul 25, 2017, at 5:03 PM, Michael Sumner notifications@github.com wrote:

@cpsievert I understand the reprojection problem with geom_raster, I'm saying the way to do that best is to use the facilities designed for raster - but also apply some tricks to make the transition between array and data frame forms a little easier for developers. This is the space I generally want to see more exploration in, and as ever I'm happy to do that with anyone.

https://gist.github.com/mdsumner/573ec70014df177baa2d1df7e55e1943

plotdap is already using projectRaster so I the extra conversion to polygons is not required. If you are interested I'll craft a PR for plotdap along those lines. The resampling applied by projectRaster is by bilinear interpolation, and it may be better to control that more explicitly to allow aggregate and other methods instead, similar to the piped example above but with an extra inverse coord transformation.

@rmendels I strongly sympathize, it's an absolute maze this stuff - at the moment it's really best to learn the available componentry rather than rely on these higher level tools. There's no language for getting between ggplot2 and spatial forms, but there are very efficient and high fidelity ways to go about it.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

---

"The contents of this message do not reflect any position of the U.S. Government or NOAA."

---

Roy Mendelssohn Supervisory Operations Research Analyst NOAA/NMFS Environmental Research Division Southwest Fisheries Science Center Note new street address 110 McAllister Way Santa Cruz, CA 95060 Phone: (831)-420-3666 Fax: (831) 420-3980 e-mail: Roy.Mendelssohn@noaa.gov www: http://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill." "From those who have been given much, much will be expected" "the arc of the moral universe is long, but it bends toward justice" -MLK Jr.

[mdsumner]

@rmendels it's missing data, easily filtered out - I definitely intend to explore the platform-specific theory of what device on what OS for the actual problem you saw. It's just not related to plotdap really and probably not related to sf is my hunch. But getting on a Mac is not easy for me. I suggest we split this and focus on getting the plots you want elsewhere. Still, it's helpful for anyone who wanted to follow this trail to be able to, so I personally don't really mind where it happens

[LDalby]

Just did a test on windows of the gist posted by @mdsumner. I see no striping in the resulting plot.

-Lars

``` Session info --------------------------------------------------------------------------- setting value version R version 3.4.1 (2017-06-30) system x86_64, mingw32 ui RStudio (1.1.299) language (EN) collate Danish_Denmark.1252 tz Europe/Paris date 2017-07-26 Packages ------------------------------------------------------------------------------- package * version date source assertthat 0.2.0 2017-04-11 CRAN (R 3.4.0) base * 3.4.1 2017-06-30 local bindr 0.1 2016-11-13 CRAN (R 3.4.1) bindrcpp 0.2 2017-06-17 CRAN (R 3.4.1) broom 0.4.2 2017-02-13 CRAN (R 3.4.0) cellranger 1.1.0 2016-07-27 CRAN (R 3.4.0) colorspace 1.3-2 2016-12-14 CRAN (R 3.4.0) compiler 3.4.1 2017-06-30 local curl 2.7 2017-06-26 CRAN (R 3.4.1) data.table 1.10.4 2017-02-01 CRAN (R 3.4.0) datasets * 3.4.1 2017-06-30 local DBI 0.7 2017-06-18 CRAN (R 3.4.1) devtools 1.13.2 2017-06-02 CRAN (R 3.4.1) digest 0.6.12 2017-01-27 CRAN (R 3.4.0) dplyr * 0.7.1 2017-06-22 CRAN (R 3.4.1) forcats 0.2.0 2017-01-23 CRAN (R 3.4.0) foreign 0.8-69 2017-06-22 CRAN (R 3.4.1) ggplot2 * 2.2.1.9000 2017-07-26 Github (tidyverse/ggplot2@1477187) git2r 0.18.0 2017-01-01 CRAN (R 3.4.0) glue 1.1.1 2017-06-21 CRAN (R 3.4.1) graphics * 3.4.1 2017-06-30 local grDevices * 3.4.1 2017-06-30 local grid 3.4.1 2017-06-30 local gtable 0.2.0 2016-02-26 CRAN (R 3.4.0) haven 1.1.0 2017-07-09 CRAN (R 3.4.1) hms 0.3 2016-11-22 CRAN (R 3.4.0) hoardr 0.2.0 2017-05-10 CRAN (R 3.4.1) httr 1.2.1 2016-07-03 CRAN (R 3.4.0) jsonlite 1.5 2017-06-01 CRAN (R 3.4.1) knitr 1.16 2017-05-18 CRAN (R 3.4.0) labeling 0.3 2014-08-23 CRAN (R 3.4.0) lattice 0.20-35 2017-03-25 CRAN (R 3.4.1) lazyeval 0.2.0 2016-06-12 CRAN (R 3.4.0) lubridate 1.6.0 2016-09-13 CRAN (R 3.4.0) magrittr 1.5 2014-11-22 CRAN (R 3.4.0) maps 3.2.0 2017-06-08 CRAN (R 3.4.1) memoise 1.1.0 2017-04-21 CRAN (R 3.4.0) methods * 3.4.1 2017-06-30 local mnormt 1.5-5 2016-10-15 CRAN (R 3.4.0) modelr 0.1.0 2016-08-31 CRAN (R 3.4.0) munsell 0.4.3 2016-02-13 CRAN (R 3.4.0) ncdf4 1.16 2017-04-01 CRAN (R 3.4.0) nlme 3.1-131 2017-02-06 CRAN (R 3.4.1) parallel 3.4.1 2017-06-30 local pkgconfig 2.0.1 2017-03-21 CRAN (R 3.4.1) plotdap * 0.0.1 2017-07-26 Github (ropensci/plotdap@3f167c1) plyr 1.8.4 2016-06-08 CRAN (R 3.4.0) psych 1.7.5 2017-05-03 CRAN (R 3.4.0) purrr * 0.2.2.2 2017-05-11 CRAN (R 3.4.0) R6 2.2.2 2017-06-17 CRAN (R 3.4.1) rappdirs 0.3.1 2016-03-28 CRAN (R 3.4.1) raster * 2.5-8 2016-06-02 CRAN (R 3.4.0) Rcpp 0.12.11 2017-05-22 CRAN (R 3.4.0) readr * 1.1.1 2017-05-16 CRAN (R 3.4.0) readxl 1.0.0 2017-04-18 CRAN (R 3.4.0) rerddap * 0.4.2 2017-05-12 CRAN (R 3.4.1) reshape2 1.4.2 2016-10-22 CRAN (R 3.4.0) rgdal 1.2-8 2017-07-01 CRAN (R 3.4.1) rlang 0.1.1 2017-05-18 CRAN (R 3.4.0) rvest 0.3.2 2016-06-17 CRAN (R 3.4.0) scales 0.4.1.9002 2017-07-06 Github (hadley/scales@6db7b6f) sf 0.5-2 2017-07-12 CRAN (R 3.4.1) sp * 1.2-5 2017-06-29 CRAN (R 3.4.1) stats * 3.4.1 2017-06-30 local stringi 1.1.5 2017-04-07 CRAN (R 3.4.0) stringr 1.2.0 2017-02-18 CRAN (R 3.4.0) tibble * 1.3.3 2017-05-28 CRAN (R 3.4.1) tidyr * 0.6.3 2017-05-15 CRAN (R 3.4.0) tidyverse * 1.1.1 2017-01-27 CRAN (R 3.4.1) tools 3.4.1 2017-06-30 local udunits2 0.13 2016-11-17 CRAN (R 3.4.0) units 0.4-5 2017-06-15 CRAN (R 3.4.1) utils * 3.4.1 2017-06-30 local withr 1.0.2 2016-06-20 CRAN (R 3.4.0) xml2 1.1.1 2017-01-24 CRAN (R 3.4.0) yaml 2.1.14 2016-11-12 CRAN (R 3.4.0) ```

[edzer]

Thanks @LDalby , that was to be expected, the question was rather whether windows shows striping on the first example of this issue, where geom_raster is not used:

library(rerddap)
library(plotdap)
library(ggplot2)
sstInfo <- info('jplMURSST41')
murSST <- griddap(sstInfo, latitude = c(22., 51.), longitude = c(-140., -105), time = c('last','last'), fields = 'analysed_sst')
add_griddap(plotdap(), murSST, ~analysed_sst, maxpixels = 50000)

It takes a while to plot; pls show us the output it gives you.

[LDalby]

Ahh, okay :-)

Those lines produces this:

[edzer]

Interesting - thanks! @hadley does any of this look familiar to you? Summary of the problem: plotting raster cells as small polygons (because they are no longer exactly a raster after reprojection) cause fine striping on mac (first post) and windows (the one above), but not on linux.

[tim-salabim]

To throw even more in the mix, here's what I get with @cpsievert 's meuse example on Windows

> sessionInfo()
R version 3.4.0 Patched (2017-05-19 r72718)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_2.2.1.9000 sf_0.5-2           raster_2.5-8       sp_1.2-5          

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12      lattice_0.20-35   digest_0.6.12     grid_3.4.0        plyr_1.8.4        gtable_0.2.0      DBI_0.6-1         magrittr_1.5     
 [9] units_0.4-4       scales_0.4.1.9002 rlang_0.1.1       lazyeval_0.2.0    labeling_0.3      rgdal_1.2-8       tools_3.4.0       udunits2_0.13    

[hadley]

IIRC this is a graphics device, and I think you can work around it by setting the (border) colourof the tiles.

[edzer]

If by "you" you mean sf, then I wonder how: e.g. here we see that sf::st_as_grob methods only take care of where to draw, not how (which colours, borders or not).

[mdsumner]

I don't have the capacity today to ensure all aspects of the environments are the same, but I believe this is enough to show that it's clearly a platform-specific effect for the edges of polygons.

Essentially, the platform dictates the difference between the default colour and the effect seen with a literal colour = NA. On Mac it's visible as edges on the polygons.

Mac: http://rpubs.com/cyclemumner/294524

Linux: http://rpubs.com/cyclemumner/294523

Windows: http://rpubs.com/cyclemumner/294528

If this is not the cause of the effect reported by @rmendels then I welcome clarification, and I apologize for the noise.

This (I think) is symptomatic of a much bigger gap in our spatial-graphics capability, and it's very topical for me. Feel free to ignore the following, but I welcome any discussion on the broader topic outside this context.

(I maintain that this a "bad" way to store and work with raster data, and if anyone is interested I'm working hard at finding better ways to do this, and I'm very happy to help implement improvements where they can help. I already know how to do it well, but wrapping it up in user-front-end is the main challenge and is what I most require assistance with. It requires grammars involving more than one table, or the use of hidden attributes used "as late as possible" to generate the mostly redundant required coordinates. With that gap in our facility, it's still better to use the single-table form with geom_raster ** as while that is limiting and wasteful, redundancy-wise it's far more efficient than converting to polygons. Coordinate transformations are supported by raster, and this is already being done in plotdap so I think that's what should be used. This is a big problem and it needs a multi-headed and planned focus put on it, as this complicated thread above illustrates.).

** (the single-centre coordinate form with auto-inference of the grain and offset, this gets interpreted and "turned back into a matrix" for rasterGrob)

[edzer]

Great examples! In the Mac one, the ugly (thin, white) striping obtained by

ggplot(p_sf, aes(fill = value)) + geom_sf(colour=NA)

disappears with

ggplot(p_sf, aes(fill = value, colour = value)) + geom_sf()

I suggest to move the "how to plot rasters and projected rasters with R" discussion to another area, for instance https://github.com/r-spatial/discuss or https://github.com/r-spatial/stars . Maybe we can come up with good guidelines/best practices, but we haven't seen all the use cases yet, and we shouldn't try here.

[edzer]

... which led to this patch in plotdap.

[trotsiuk]

The issue seems to persist on Mac if one add alpha = ...

library(raster)
library(sf)
library(ggplot2)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- st_as_sf(rasterToPolygons(r))
ggplot() + geom_sf(data = s, aes(fill = test, colour = test), alpha = 0.8)

Session Info

``` R version 3.4.3 (2017-11-30) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS High Sierra 10.13.6 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] ggplot2_3.0.0 sf_0.6-3 raster_2.6-7 sp_1.3-1 loaded via a namespace (and not attached): [1] Rcpp_0.12.18 pillar_1.3.0 compiler_3.4.3 plyr_1.8.4 bindr_0.1.1 [6] class_7.3-14 tools_3.4.3 digest_0.6.15 tibble_1.4.2 gtable_0.2.0 [11] lattice_0.20-35 pkgconfig_2.0.1 rlang_0.2.1.9000 DBI_1.0.0 rstudioapi_0.7 [16] yaml_2.2.0 rgdal_1.2-18 bindrcpp_0.2.2 spData_0.2.9.0 e1071_1.6-8 [21] withr_2.1.2 dplyr_0.7.6 classInt_0.2-3 grid_3.4.3 tidyselect_0.2.4 [26] glue_1.3.0 R6_2.2.2 purrr_0.2.5 udunits2_0.13 magrittr_1.5 [31] scales_1.0.0.9000 units_0.5-1 assertthat_0.2.0 colorspace_1.3-2 labeling_0.3 [36] lazyeval_0.2.1 munsell_0.5.0 crayon_1.3.4 ```

[edzer]

I see that too on ubuntu.

[trotsiuk]

I think it might be coming from geom_polygon() as there is the same issue. alpha only apply on fill and not on color for the geom_polygon (also) and I think this is how it works with st_geom.

[ComputationalEpi]

Came across this thread when having similar issues and was frustrated that there was no solution. I found if you use CairoPDF() from the Cairo package that this problem is eliminated.