walkerke / mapgl

R interface to Mapbox GL JS v3 and Maplibre GL JS
https://walker-data.com/mapgl
Other
91 stars 5 forks source link

`add_image_source()` doesn't accept colour tab in `spatRaster` #21

Open RWParsons opened 4 months ago

RWParsons commented 4 months ago

There are a couple problem that I've identified in add_image_source()

Here is an example using a sample of data that I'm working with:

test_data <- structure(list(x = c(146.704635285663, 143.470715055695, 145.815307222422, 
                                  147.917355371901, 145.060725835429, 144.117499101689, 148.186848724398, 
                                  139.805605461732, 148.456342076895, 149.04922745239, 140.641034854474, 
                                  144.575637800934, 147.351419331656, 143.335968379447, 143.335968379447, 
                                  141.314768235717, 145.545813869925, 152.714337046353, 147.593963348904, 
                                  150.747035573123, 146.327344592167, 148.779734099892, 148.96837944664, 
                                  147.432267337406, 144.198347107438, 150.908731584621, 144.198347107438, 
                                  142.931728350701, 144.548688465685, 144.710384477183, 148.213798059648, 
                                  141.422565576716, 147.782608695652, 140.506288178225, 143.01257635645, 
                                  140.425440172476, 143.282069708947, 144.791232482932, 140.398490837226, 
                                  151.636363636364, 151.98670499461, 144.494789795185, 141.98850161696, 
                                  144.872080488681, 143.01257635645, 142.716133668703, 145.680560546173, 
                                  144.602587136184, 144.683435141933, 142.177146963708), y = c(-22.270970176069, 
                                                                                               -21.6241861300755, -22.7560582105641, -22.567412863816, -28.6040639597557, 
                                                                                               -20.3306180380884, -20.3306180380884, -18.6597592526051, -23.5375889328063, 
                                                                                               -25.1545490477902, -20.573162055336, -15.9648257276321, -20.573162055336, 
                                                                                               -18.0399245418613, -16.6385591088753, -21.0043514193317, -24.6964103485447, 
                                                                                               -26.6098131512756, -28.3076212720086, -25.0467517067912, -21.5433381243263, 
                                                                                               -23.3758929213079, -23.5375889328063, -27.1218505210205, -17.1775458138699, 
                                                                                               -25.3162450592885, -21.1121487603306, -23.1063995688106, -23.1872475745598, 
                                                                                               -16.8541537908732, -20.4384153790873, -23.4567409270571, -24.9120050305426, 
                                                                                               -21.1929967660798, -20.8426554078333, -24.2382716492993, -16.9080524613726, 
                                                                                               -27.2296478620194, -22.7021595400647, -24.8581063600431, -26.5559144807761, 
                                                                                               -26.4481171397772, -23.1333489040604, -23.4028422565577, -14.267017606899, 
                                                                                               -18.471113905857, -18.4441645706073, -24.0765756378009, -19.2795939633489, 
                                                                                               -19.6568846568451), pred = c(305.917283580099, 357.147308791887, 
                                                                                                                            376.899080177541, 311.293558669515, 385.302194465796, 311.406970170659, 
                                                                                                                            213.397097913667, 381.62892998555, 380.620131247899, 323.160634223337, 
                                                                                                                            353.569875088294, 252.367408752526, 200.205388159247, 329.665758425744, 
                                                                                                                            303.616701525762, 347.401323384892, 388.910507958226, 163.814361430662, 
                                                                                                                            361.271026317014, 254.712761913955, 279.277148678091, 372.939247974351, 
                                                                                                                            369.336082811912, 345.740485904605, 250.928162986335, 246.393765060807, 
                                                                                                                            337.069420587925, 380.044721245275, 407.590072013831, 221.797793265574, 
                                                                                                                            222.250923169777, 362.431972284566, 380.155056796395, 368.03688613144, 
                                                                                                                            340.701795340508, 395.745150086577, 307.098348063308, 377.360232941649, 
                                                                                                                            368.020435732177, 320.938285611497, 174.648842951226, 374.875642089735, 
                                                                                                                            362.421003870795, 406.694805719537, 385.248720216423, 349.138641118029, 
                                                                                                                            178.077076149705, 400.417789912908, 248.847855237357, 350.543773979372
                                                                                               ), id = c(119188L, 129383L, 110793L, 114199L, 4932L, 149130L, 
                                                                                                         149281L, 170610L, 97078L, 66001L, 145388L, 194585L, 145637L, 
                                                                                                         177759L, 190621L, 138950L, 75096L, 38399L, 9545L, 68258L, 130767L, 
                                                                                                         99973L, 97097L, 29644L, 186468L, 62739L, 137392L, 104512L, 103145L, 
                                                                                                         189183L, 147669L, 98269L, 70863L, 136007L, 141475L, 83781L, 188708L, 
                                                                                                         27749L, 111544L, 72086L, 39270L, 40788L, 104003L, 99352L, 202832L, 
                                                                                                         172885L, 173304L, 87022L, 163546L, 158565L), col = c("#FA4A29", 
                                                                                                                                                              "#E41E1D", "#E1191C", "#F74627", "#E1191D", "#F74627", "#FE9E43", 
                                                                                                                                                              "#E1191D", "#E1191D", "#F23D24", "#E6221D", "#FD8238", "#FEA647", 
                                                                                                                                                              "#F03823", "#FA4C29", "#E8291F", "#E0191D", "#FEBD57", "#E31A1C", 
                                                                                                                                                              "#FD7F37", "#FD6730", "#E2191C", "#E21A1C", "#E92A1F", "#FD8339", 
                                                                                                                                                              "#FD873A", "#ED3221", "#E1191D", "#DF181D", "#FE9941", "#FE9841", 
                                                                                                                                                              "#E31A1C", "#E1191D", "#E21A1C", "#EB2F20", "#E0191D", "#F94928", 
                                                                                                                                                              "#E1191C", "#E21A1C", "#F33F25", "#FEB650", "#E2191C", "#E31A1C", 
                                                                                                                                                              "#DF181D", "#E1191D", "#E8271E", "#FEB34D", "#DF181D", "#FD8539", 
                                                                                                                                                              "#E7251E")), row.names = c(NA, -50L), class = c("tbl_df", "tbl", 
                                                                                                                                                                                                              "data.frame"))

tb <- data.frame(value = test_data$id, col = test_data$col)

test_raster <- tidyterra::as_spatraster(dplyr::select(test_data, x, y, id), crs = 4326)

terra::coltab(test_raster) <- tb

plot(test_raster) # works

maplibre() |>
  # add rasters
  add_image_source(
    id = "acute_raster_data",
    data = test_raster
  ) |>
  add_raster_layer(
    id = "acute_raster_layer",
    source = "acute_raster_data"
  )

This is just a sample of 50 rows from the input data. The full data is of an Australian state and looks fine when running plot(test_raster):

image

One of the issues is that this raster has colour layer that add_image_source() should accept. However, after it is projected to EPSG:4326 here: https://github.com/walkerke/mapgl/blob/7f3e2a5932dd466f7b5e4874127cbb40d99580d3/R/sources.R#L161

This causes it to lose the it's attribute and terra::has.colors(data) then returns FALSE. This could arguably being a feature request to terra to but I think it could be solved if the coltb is created before the data are projected in add_image_source().

The second problem I found is that even if I skip the projection, with these data, I get an error when running: https://github.com/walkerke/mapgl/blob/7f3e2a5932dd466f7b5e4874127cbb40d99580d3/R/sources.R#L190

I don't think the rgb values are out of range but it fails to write the image and then nothing is shown on the map.

If you would like, I can make a PR to solve the first issue with the projection but I'm not sure what to do about the problem with writeRaster (and you're probably much more familiar with this).

Cheers, Rex

RWParsons commented 4 months ago

update: I was able to get the writeRaster step to work by making tb a maximum of 255 rows (only 255 possible values for the outcome variable in the spatRaster). This still gets messed up by the projection, though. Happy to push a fix if you'd like, @walkerke.

walkerke commented 4 months ago

Yes a PR would be great!

walkerke commented 4 months ago

In debugging this, I'm also noticing that terra::project() resamples the raster, so we have to project before doing the data rescaling. I think I'm getting closer to an implementation that'll work here.

RWParsons commented 4 months ago

Right - I've been struggling to get it working to do the projection after and reliably add the colour tab back in. Not quite sure what the best approach will be to address this issue.

walkerke commented 4 months ago

I modified the code to throw a warning for now that the existing color table won't be used. Because you no longer have a 1:1 correspondence between cells after projection, we'd need to establish the old color table, project, then associate the value in the new raster with the corresponding nearest value in the old color table. I haven't quite gotten this working yet.