Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
rmendels opened this issue Jul 24, 2017 · 39 comments
Closed

Comments

@rmendels
Copy link

rmendels commented Jul 24, 2017

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

plotdapgg
plotdapbase

@edzer
Copy link
Member

edzer commented Jul 24, 2017

Thanks! Strange enough, I see
x

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
Copy link
Member

edzer commented Jul 24, 2017

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

@rmendels
Copy link
Author

rmendels commented Jul 24, 2017 via email

@edzer
Copy link
Member

edzer commented Jul 24, 2017

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

@edzer edzer closed this as completed Jul 24, 2017
@edzer
Copy link
Member

edzer commented Jul 24, 2017

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

@cpsievert
Copy link
Contributor

cpsievert commented Jul 24, 2017

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

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
Copy link
Author

rmendels commented Jul 24, 2017 via email

@rmendels
Copy link
Author

rmendels commented Jul 24, 2017 via email

@mdsumner
Copy link
Member

mdsumner commented Jul 25, 2017

@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
Copy link
Author

rmendels commented Jul 25, 2017 via email

@edzer
Copy link
Member

edzer commented Jul 25, 2017

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
Copy link
Member

mdsumner commented Jul 25, 2017

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
Copy link
Author

rmendels commented Jul 25, 2017 via email

@rmendels
Copy link
Author

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:

rplot

Michael Sumner asked what happens if I use the following:

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

I get:

rplot

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

rplot

@cpsievert
Copy link
Contributor

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

image

@rmendels
Copy link
Author

rmendels commented Jul 25, 2017 via email

@tim-salabim
Copy link
Member

tim-salabim commented Jul 25, 2017

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

screenshot at 2017-07-25 17 20 15

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
Copy link
Author

rmendels commented Jul 25, 2017 via email

@mdsumner
Copy link
Member

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
Copy link
Contributor

cpsievert commented Jul 25, 2017

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
Copy link
Author

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
Copy link
Member

mdsumner commented Jul 26, 2017

@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
Copy link
Author

rmendels commented Jul 26, 2017 via email

@rmendels
Copy link
Author

rmendels commented Jul 26, 2017 via email

@mdsumner
Copy link
Member

@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
Copy link

LDalby commented Jul 26, 2017

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

-Lars

image

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
Copy link
Member

edzer commented Jul 26, 2017

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
Copy link

LDalby commented Jul 26, 2017

Ahh, okay :-)

Those lines produces this:
image

@edzer
Copy link
Member

edzer commented Jul 26, 2017

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
Copy link
Member

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

rplot

> 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
Copy link
Contributor

hadley commented Jul 26, 2017

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

@edzer
Copy link
Member

edzer commented Jul 26, 2017

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
Copy link
Member

mdsumner commented Jul 27, 2017

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
Copy link
Member

edzer commented Jul 27, 2017

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
Copy link
Member

edzer commented Jul 28, 2017

... which led to this patch in plotdap.

karawoo pushed a commit to tidyverse/ggplot2 that referenced this issue Jul 28, 2017
* addresses #2119

* tidy graticule fixes

* ggplot2 side of issue r-spatial/sf#396

* tabs -> spaces

* add ndiscr to docs

* fix #2200

* attempt to fix #2060

All cases were in sf.R a geometry column is address with x$geometry, ggplot2 made the wrong assumption that the geometry column has a fixed name. I replaced this in certain instances, where the data are already pretty transformed and no longer have properties of sf objects, with a fixed position, i.e. x[[1]], which seems to work.

* fixes r-spatial/sf#438

* address review comments

* fix break on geom_raster, objects without list-column

see https://gist.github.com/mdsumner/573ec70014df177baa2d1df7e55e1943 for the case that this PR fixes

* tidy up

* trying @karawoo's suggestion

* adds some sf tests

* tidy further
@trotsiuk
Copy link

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)

rplot

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
Copy link
Member

edzer commented Aug 24, 2018

I see that too on ubuntu.

@trotsiuk
Copy link

trotsiuk commented Aug 24, 2018

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.

elizabethng added a commit to elizabethng/predator-diets that referenced this issue Dec 8, 2019
@ComputationalEpi
Copy link

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

9 participants