search

rsbivand / eqc23_talk

Talk at ECTQG2023, Braga ( Portugal) 14 -17 September 2023
Other
0 stars 0 forks source link

RFC talk draft #1

Open rsbivand opened 1 year ago

rsbivand commented 1 year ago

@rCarto @HuguesPecout @geocaruso @dieghernan

I've put a very preliminary presentation draft up on https://rsbivand.github.io/eqc23_talk/, repo https://github.com/rsbivand/eqc23_talk (here). I'd be very grateful for comments.

The presentation is only a first approach, but I'd like to develop this towards a vignette, and would also like to bring in package maintainers of packages using class intervals to see how to provide support to users.

geocaruso commented 1 year ago

Hi Roger, Great idea. Will have some further look at get back to you.

I am not sure how useful but I have a repo for mapping Luxembourg population per municipality or using the EU ref grid (1km). I believe it can bring 2 discussion points:

Besides you’ll see I am basically using ggplot for mapping with a function I have named ggplot.themap() (see R folder and its application to the Lux case in ggplot.themap.lux). My idea was originally to use the dots to pass parameters to classInt but I retracted because of some bugs I could not understand with the “fixed” method.

Some output with varying discretisation methods here: https://github.com/geocaruso/cartolux/blob/main/output/Lux_map_grid_pop.pdf

In general it may also just show how a rather basic user (not a proper developer) makes use of classint in cartographic context with public data only.

Best Geoffrey

On 17 Aug 2023, at 13:46, Roger Bivand @.***> wrote:

@rCarto https://github.com/rCarto @HuguesPecout https://github.com/HuguesPecout @geocaruso https://github.com/geocaruso @dieghernan https://github.com/dieghernan I've put a very preliminary presentation draft up on https://rsbivand.github.io/eqc23_talk/, repo https://github.com/rsbivand/eqc23_talk (here). I'd be very grateful for comments.

The presentation is only a first approach, but I'd like to develop this towards a vignette, and would also like to bring in package maintainers of packages using class intervals to see how to provide support to users.

— Reply to this email directly, view it on GitHub https://github.com/rsbivand/eqc23_talk/issues/1, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASTDGZRUSPUKQOX6QXLUOETXVYAAHANCNFSM6AAAAAA3T6FBXY. You are receiving this because you were mentioned.

rsbivand commented 1 year ago

@geocaruso Thanks, yes, plenty of useful points to consider. I'll update the draft closer to the meeting.

geocaruso commented 1 year ago

As I needed to do something similar today, I attempted to use classInt for a spatraster with ggplot2 (following the line of what I do with sf's). Well I believe it is quite a painful experience. tidyterra helps a lot with continuous values mapping but as far as I could see does not work with discrete values

Here my code using 1km pop for Luxembourg.

`

1. Data

Grid1km_LAU2_Pop2021EU<-read.csv("https://github.com/geocaruso/cartolux/raw/main/data/Grid1km_LAU2_Pop2021EU.csv")

GrdPop<-Grid1km_LAU2_Pop2021EU[,c(2,6)] #extract cell code and pop only

2. Make raster from ref grid:

extract lower left corner, then factor 1000 for meters, and add 500m to get to centre

x_centre<-as.numeric(substr(GrdPop$CELLCODE,5,8))1000+500
y_centre<-as.numeric(substr(GrdPop$CELLCODE,10,13))1000+500

spatraster from xyz

xyz0<-cbind(x_centre,y_centre,GrdPop$Pop2021EU) xyz<-xyz0[xyz0[,3]>0,] r<-terra::rast(xyz, type="xyz") terra::crs(r) <- "epsg:3035"

3. Map pop SpatRaster with ggplot and continuous scale

EASY with tidyterra!

library(ggplot2) p<-ggplot()+ tidyterra::geom_spatraster(data = r, aes(fill = X))+ scale_fill_viridis_c(option="magma",na.value = "white",direction = -1) p

4. Map SpatRaster with ggplot and discretized scale??

NOT EASY!

cl.intvl<-classInt::classIntervals(xyz[,3], n=5, style="fisher") n.classes<-length(cl.intvl$brks)-1 cl.value<-factor(classInt::findCols(cl.intvl))

xyz2_df<-as.data.frame(cbind(xyz[,1:2],cl.value))

p2<-ggplot()+ geom_raster(data = xyz2_df, aes(x=x_centre, y=y_centre,fill = factor(cl.value)))+ tidyterra::geom_spatraster(data = r, aes(fill = factor(cl.value)))+ scale_fill_brewer(palette="Spectral", direction=-1) p2

Compare p and p2

Conclusion: Quite ugly ;-)

rsbivand commented 1 year ago

I think we'd need to look at terra::plot as well as support in tmap3 and forthcoming tmap4, maybe others. I see that mapsf has some support too: https://riatelab.github.io/mapsf/reference/mf_raster.html. I stay away from ggplot2 unless forced; the analysis in How does ggplot do this? in https://rsbivand.github.io/BAN422_V20/ban422_v20_tues.html, which takes functions calling functions to extremities IMO.

geocaruso commented 1 year ago

I think you are right about tmap vs gglot. I can’t really remember why I quit tmap, especially since the leaflet mode is very practical.

As to the SpatRaster I think the logical way to go is to feed the rcl argument of terra::classify with a from to and new value matrix or with the breaks from the classInt object

Geoffrey

Le mar. 22 août 2023 à 21:19, Roger Bivand @.***> a écrit :

I think we'd need to look at terra::plot as well as support in tmap3 and forthcoming tmap4, maybe others. I see that mapsf has some support too: https://riatelab.github.io/mapsf/reference/mf_raster.html. I stay away from ggplot2 unless forced; the analysis in How does ggplot do this? in https://rsbivand.github.io/BAN422_V20/ban422_v20_tues.html, which takes functions calling functions to extremities IMO.

— Reply to this email directly, view it on GitHub https://github.com/rsbivand/eqc23_talk/issues/1#issuecomment-1688788879, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASTDGZV5RAFE4Q52KKMOMC3XWUA5TANCNFSM6AAAAAA3T6FBXY . You are receiving this because you were mentioned.Message ID: @.***>

dieghernan commented 8 months ago

Very late on this, as contributor to classInt didn't have any comment but as developer of tidyterra...

As I needed to do something similar today, I attempted to use classInt for a spatraster with ggplot2 (following the line of what I do with sf's). Well I believe it is quite a painful experience. tidyterra helps a lot with continuous values mapping but as far as I could see does not work with discrete values @geocaruso

I think plotting SpatRasters with discrete values is also very straighforward, see:

Grid1km_LAU2_Pop2021EU <- read.csv("https://github.com/geocaruso/cartolux/raw/main/data/Grid1km_LAU2_Pop2021EU.csv")

GrdPop <- Grid1km_LAU2_Pop2021EU[, c(2, 6)] # extract cell code and pop only

# extract lower left corner, then factor 1000 for meters, and add 500m to get to centre
x_centre <- as.numeric(substr(GrdPop$CELLCODE, 5, 8)) * 1000 + 500
y_centre <- as.numeric(substr(GrdPop$CELLCODE, 10, 13)) * 1000 + 500

xyz0 <- cbind(x_centre, y_centre, GrdPop$Pop2021EU)
xyz <- xyz0[xyz0[, 3] > 0, ]
r <- terra::rast(xyz, type = "xyz")
terra::crs(r) <- "epsg:3035"

# 3. Map pop SpatRaster with ggplot and continuous scale
# EASY with tidyterra!

library(ggplot2)

p <- ggplot() +
  tidyterra::geom_spatraster(data = r) +
  scale_fill_viridis_c(option = "magma", na.value = "white", direction = -1)
p

For making a SpatRaster discrete, use terra::classify() with the breaks of classInt::classIntervals() so you can create a categorical layer.


# 4. Map SpatRaster with ggplot and discretized scale
cl.intvl <- classInt::classIntervals(xyz[, 3], n = 5, style = "fisher")

# Classify the SpatRaster into factors
r$classified <- terra::classify(r, cl.intvl$brks,
  include.lowest = TRUE,
  brackets = TRUE
)

When we are done, just include in aes(fill=) the name of the new layer with the factors.

# Plot SpatRaster with ggplot2
p2 <- ggplot() +
  tidyterra::geom_spatraster(data = r, aes(fill = classified)) + # Note the layer we use
  scale_fill_brewer(palette = "Spectral", direction = -1, na.translate = FALSE) + # na.translate remove NA for legend
  labs(fill = "Using factors")

p2

If we prefer numeric factors...


# New layer
integs <- terra::values(r$classified)
r$num_as_fact <- as.factor(as.integer(integs) + 1)

p3 <- ggplot() +
  tidyterra::geom_spatraster(data = r, aes(fill = num_as_fact)) +
  scale_fill_brewer(palette = "Spectral", direction = -1, na.translate = FALSE) +
  labs(fill = "A number as factor")

p3

Created on 2024-02-26 with reprex v2.1.0

rsbivand commented 8 months ago

@dieghernan thanks!

geocaruso commented 7 months ago

Thanks @dieghernan. I see now I overlooked terra::classify which can make us the classint breaks quickly indeed. And through a factoring of terra::values for integers rather than breaks. Nice! Much better than getting back to xyz.