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

implement new df to pts to raster with possibly different xres() and yres() #3

Open
bbest opened this issue Jul 19, 2023 · 0 comments

Comments

@bbest
Copy link
Contributor

bbest commented Jul 19, 2023

TODO: NEW

Per noaa-onms/climate-dashboard: /scripts/get_data.R#L182-L233

        # get geospatial attributes
        a <- ed$alldata$NC_GLOBAL |>
          filter(
            attribute_name |> str_starts("geospatial_"),
            data_type == "double") |>
          select(attribute_name, value)
        g <- setNames(as.numeric(a$value), a$attribute_name) |> as.list()
        lon_half <- g$geospatial_lon_resolution/2
        lat_half <- g$geospatial_lat_resolution/2

        # setup raster with potentially different xres() and yres()
        r_template <- rast(
          xmin       = min(x$lon) - lon_half,
          xmax       = max(x$lon) + lon_half,
          ymin       = min(x$lat) - lat_half,
          ymax       = max(x$lat) + lat_half,
          resolution = c(
            g$geospatial_lon_resolution,
            g$geospatial_lat_resolution),
          crs = "epsg:4326")

        df_to_rast <- function(df, r_template){
          # data frame to points
          p <- df |>
            select(lon, lat, all_of(ed_row$var)) |>
            st_as_sf(
              coords = c("lon", "lat"),
              crs    = 4326)
          # points to raster
          rasterize(p, r_template, field = ed_row$var)
        }

        x <- x |>
          group_by(date) |>
          nest(data = all_of(c("lon", "lat", ed_row$var))) |>
          mutate(
            r = map(data, df_to_rast, r_template))

        stk <- rast(x$r)
        names(stk) <- glue("{ed_row$var}_{as.character(x$date) |> str_replace_all('-','.')}")
        crs(stk) <- "EPSG:4326"

CURRENT

extractr/R/read.R

Lines 334 to 390 in 3b8feb6

if (all(c("lon", "lat") %in% colnames(nc$data))){
d <- tibble(nc$data) %>%
mutate(
# round b/c of uneven intervals
# unique(tbl$lon) %>% sort() %>% diff() %>% unique() %>% as.character()
# 0.0499954223632812 0.0500030517578125
# TODO: inform Maria/Joaquin about uneven intervals
lon = round(lon, 4),
lat = round(lat, 4),
date = as.Date(time, "%Y-%m-%dT12:00:00Z"))
} else if (all(c("longitude", "latitude") %in% colnames(nc$data))){
d <- tibble(nc$data) %>%
mutate(
lon = round(longitude, 4),
lat = round(latitude, 4),
# lon = longitude,
# lat = latitude,
date = as.Date(time, "%Y-%m-%dT12:00:00Z"))
} else {
stop("Expected lon/lat or longitude/latitude in ERDDAP dataset.")
}
d_sp <- d |>
select(lon, lat, !!ed_var)
sp::coordinates(d_sp) <- ~ lon + lat
# x0 <- x
# sp::gridded(x) <- T
# r <- raster::raster(x, layer = ed_var)
if (ed_var == "chlorophyll"){
g <- sp::points2grid(d_sp, tolerance = 0.0243902)
} else {
g <- sp::points2grid(d_sp)
}
# g <- sp::points2grid(x, tolerance = 1e-05)
# g <- sp::points2grid(x), tolerance = 1e-05)
# cx <- range(diff(sort(unique(d$lon))))
# # dx <- range(diff(cx))
# cy <- range(diff(sort(unique(d$lat))))
# # dy <- range(diff(cy))
# diff(c(unique(cx), unique(cy)))
#
# diff(cx), diff(cy)
# tibble(
# cx = cx,
# cy = cy) |>
# expand(cx, cy) |>
# mutate(dif = )
#
# dy <- range(diff(sort(unique(d$lat))))
# tol <- max(max(c(dx,dy)))
# g <- try(sp::points2grid(x))
# g <- sp::points2grid(x, tolerance = tol)
r <- raster::raster(g)
idx <- raster::cellFromXY(r, sp::coordinates(d_sp))
r[idx] <- d_sp[[ed_var]]
# plot(r)
raster::crs(r) <- 4326
r

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

1 participant