warp example #214
Comments
|
and, does the json info for geolocation arrays, GCPs, RPCs, reflect the actual footprint? has implications for the wider use of the footprint utility (which is geared at nearblack collars atm) |
|
can now do dsn = 'NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/LST_in.nc":LST'
geol = c('vrt://NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/geodetic_in.nc":longitude_in?unscale=true&ot=Float64',
'vrt://NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/geodetic_in.nc":latitude_in?unscale=true&ot=Float64')
vrt <- vapour::vapour_vrt(dsn, geolocation = geol)
x <- vapour::gdal_raster_data(vrt, target_dim = c(512, 0))
|
|
can do this to get through sf: sf::gdal_utils("warp", vrt, "/vsimem/warp.tif", options = c("-ts", "512", "0"))
terra::plot(terra::rast("/vsimem/warp.tif"))
not sure if there's a way to add GEOLOCATION domain metadata to a utils translate call ... |
|
Hi!
There are missing pixels. And there is also important information like quality, cloud mask, and orphan pixels. |
|
what missing pixels?? I appreciate you looking at it but there's a lot of missing context here, it's purely a note fir myself for now 🙏 this is a purely technical situation being tested, it's not intended to be aware of the meaning or intention of the data- there's things about the setting up of geolocation arrays that stars is doing internally (when it does this its way) and I want that to be transparent and easily available in GDAL too your case doesn't do the type conversion and scaling of the arrays /1e6, and doing that with the vrt:// syntax is only available in more recent GDAL (which is all I'm looking at here) |
|
I wanted to explore your approach since you referred to the stars issue I created. What I meant about the other information is that the re-gridding process needs to be carefully conducted. What if you want to do this with categorical data (as flag/quality layer) and it does simple interpolation between pixels? That is why the Sentinel 3 LST dataset is a bit problematic. Ok, since it is a note for yourself I won't bother. I understand it is about GDAL possibilities, not the data itself. The issue posted by me was referred to, so I followed:) |
|
I didn't realize it was your example, I'm sorry. Thanks for the the blue, that makes it clear. I'm not interested in the right re-gridding process for this data (yet), I want the GDAL part to work well first and not be lost in downstream R package obscurity. I was intending to come back to you with the full example, but I'm not there yet - which is why I didn't ping the original issue it was intended just for me to look at until I had it fully worked. There are things in GDAL itself that won't be fixed or (more) convenient for this until at least GDAL 3.8. I agree about the auxiliary information, but before that there is precursor stuff about getting the geoloation arrays right, and I don't like that stars hides that away - if we get it in GDAL it will work the same in all downstream languages. I'm particularly also looking at the OpenDataCube and how it would approach this, but it's early days for me there. I'm very happy to explore this with you! Apologies for being unclear and taking it away from your original post. |
|
when done at full resolution, and with the correct type and unscaling of the geolocation arrays there's no missing data (and this is purely with default nearest neighbour resampling:
Just note that syntax (the dsn = 'NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/LST_in.nc":LST'
geol = c('vrt://NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/geodetic_in.nc":longitude_in?unscale=true&ot=Float64',
'vrt://NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/geodetic_in.nc":latitude_in?unscale=true&ot=Float64')
vrt <- vapour::vapour_vrt(dsn, geolocation = geol)
vapour::vapour_raster_info(dsn)$dimension
## give no arguments and you get best-native resolution (but don't do this if the info$dim is really massive!)
x <- vapour::gdal_raster_data(vrt)
pt <- vaster::xy_from_cell(attr(x, "dimension"), attr(x, "extent"), which(is.na(x[[1]])))
points(pt, pch = 19, cex = 0.2, col = "blue")
|
|
a little bit more, I think this might work with GDAL < 3.7.0, and see I'm using /vsimem to avoid intermediary temp files - but, still I'm working in dev-GDAL and haven't tested in older versions (I'm very glad to see this works as I intend though! there's a lot more to say) also of course we should target a specific grid not just let GDAL pick a longlat extent, but see now we have the right scaling for longlat and without internal tricks in sf or stars, this is pure GDAL and would work the same in python ## to do this now with GDAL < 3.7.0 sf utils do
dsn = 'NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/LST_in.nc":LST'
geol_scaled = c('NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725/S3A_SL_2_LST.zip/geodetic_in.nc":longitude_in',
'NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725/S3A_SL_2_LST.zip/geodetic_in.nc":latitude_in')
tf_lon <- "/vsimem/lon_unscaled.tif"
tf_lat <- "/vsimem/lat_unscaled.tif"
sf::gdal_utils("translate", geol_scaled[1], tf_lon, options = c("-ot", "Float64", "-unscale"))
sf::gdal_utils("translate", geol_scaled[2], tf_lat, options = c("-ot", "Float64", "-unscale"))
vrt <- vapour::vapour_vrt(dsn, geolocation = c(tf_lon, tf_lat))
## use sf or stars or vapour, even terra::project probably
sf::gdal_utils("warp", vrt, "/vsimem/warp.tif")
x <- terra::rast("/vsimem/warp.tif")
terra::plot(x, background = "blue")## what if we do that at lower resolution (for quick look) - set one dim to 0 so GDAL figures out a good value
sf::gdal_utils("warp", vrt, "/vsimem/warp_view2.tif", options = c("-ts", c(112, 0)))
x <- terra::rast("/vsimem/warp_view2.tif")
terra::plot(x, background = "blue")Created on 2023-09-20 with reprex v2.0.2 finally one more example in a sensible map projection for the region ## to do this now with GDAL < 3.7.0 sf utils do
dsn = 'NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725//S3A_SL_2_LST.zip/LST_in.nc":LST'
geol_scaled = c('NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725/S3A_SL_2_LST.zip/geodetic_in.nc":longitude_in',
'NETCDF:"/vsizip//vsicurl/https://github.com/r-spatial/stars/files/10984725/S3A_SL_2_LST.zip/geodetic_in.nc":latitude_in')
tf_lon <- "/vsimem/lon_unscaled.tif"
tf_lat <- "/vsimem/lat_unscaled.tif"
sf::gdal_utils("translate", geol_scaled[1], tf_lon, options = c("-ot", "Float64", "-unscale"))
sf::gdal_utils("translate", geol_scaled[2], tf_lat, options = c("-ot", "Float64", "-unscale"))
vrt <- vapour::vapour_vrt(dsn, geolocation = c(tf_lon, tf_lat))
## use sf or stars or vapour, even terra::project probably
sf::gdal_utils("warp", vrt, "/vsimem/local.tif", options = c("-t_srs", "+proj=lcc +lon_0=67 +lat_0=36 +lat_1=42 +lat_2=28 +datum=WGS84"))
x <- terra::rast("/vsimem/local.tif")
x
#> class : SpatRaster
#> dimensions : 1481, 1720, 1 (nrow, ncol, nlyr)
#> resolution : 992.7009, 992.7009 (x, y)
#> extent : -825457.1, 881988.4, -786815, 683374.9 (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=lcc +lat_0=36 +lon_0=67 +lat_1=42 +lat_2=28 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#> source : local.tif
#> name : local
terra::plot(x, background = "black", col = hcl.colors(128))
m <- terra::project(do.call(cbind, maps::map(plot = F)[1:2]), to = x, from = "OGC:CRS84")
#> Warning: [project] 1972 failed transformations
lines(m)Created on 2023-09-20 with reprex v2.0.2 |
|
i'm also exploring ways to avoid that creation of VRT with vapour, it's easy enough to craft manually but it should be possible with gdal.Translate, and I'm working on that :) |
|
That looks great! Thank you for your work on this. If you need an extra pair of hands for any testing (I am not a developer, but a heavy user), count me in!:) Thanks again! |
|
cool, how are you with running docker? can run recent builds of gdal using it so we can test stuff I'm not entirely comfy with it myself yet (I build gdal from source on a machine with R already installed), but rocker is working on making it easier to fold R and gdal together |
|
I have never set up a Docker container; never needed for one. But I know it has become very popular. |
see how this looks when done as gdal intended
stars/issues/618
the problem is the lon,lat units are scaled * 1e-6
so, possibly the SRS could contain that unit conversion
The text was updated successfully, but these errors were encountered: