AustralianAntarcticDivision / raadtools

Tools for Synoptic Environmental Spatial Data
http://australianantarcticdivision.github.io/raadtools/
24 stars 3 forks source link

polarview JPEG images #126

Open mdsumner opened 3 years ago

mdsumner commented 3 years ago

Rough workflow, functions to get the jpegs/tif-tarballs, and some code to get the tif dimensions and geotransform, convert that to an extent for the JPEG

files <- polarview_files(type = "jpeg")
afile <- files$fullname[nrow(files)-1]
tarball <- raadfiles:::polarview_jpeg_tarball(afile)

#tarball <- na.omit(tarball)
## must be valid tarball paths here (not missing)
#polarview_get_geotransform(tarball)

info <- vapour::vapour_raster_info(glue::glue("/vsitar/{tarball}/{polarview_tifname(tarball)}"))

## new vapour has extent
(ex <- info$extent)
# 1757376  2063256 -1934750 -1477390

# Upper Left  ( 1757375.669,-1477389.517) (130d 3'11.07"E, 69d 5'42.72"S)
# Lower Left  ( 1757375.669,-1934749.517) (137d45' 1.58"E, 66d16'30.90"S)
# Upper Right ( 2063255.669,-1477389.517) (125d36'16.10"E, 66d56'52.03"S)
# Lower Right ( 2063255.669,-1934749.517) (133d 9'32.40"E, 64d23' 8.63"S)
# Center      ( 1910315.669,-1706069.517) (131d46' 2.94"E, 66d44'20.76"S)

## now read the JPEG and use that extent (the transform doesn't apply because jpeg and
## tiff have different resolution, same extent though it seems)

jpeg <- raster::brick(afile)
extent(jpeg) <- ex
plotRGB(jpeg)

## can also read the tarball direct
terra::rast(info$filelist[1])
class       : SpatRaster 
dimensions  : 14815, 15171, 1  (nrow, ncol, nlyr)
resolution  : 40, 40  (x, y)
extent      : 2651343, 3258183, -395393.4, 197206.6  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : S1A_EW_GRDM_1SDH_20210928T133512_33E3_S_1.tif 
name        : S1A_EW_GRDM_1SDH_20210928T133512_33E3_S_1 
mdsumner commented 2 years ago

Here's an example update with the new GDAL/vapour/helper features, it simplfies what we were doing above only a little, but is very generally applicable, and is independent of the other spatial packages in R:

library(raadfiles)
files <- polarview_files(type = "jpeg")
afile <- files$fullname[nrow(files)]
tarball <- raadfiles:::polarview_jpeg_tarball(afile)

library(whatarelief)
info1 <- vapour::vapour_raster_info(glue::glue("/vsitar/{tarball}/{raadfiles:::polarview_tifname(tarball)}"))

vrt <- vapour::vapour_vrt(afile, extent  = info1$extent, 
                          projection = info1$projection)

im <- imagery(source = vrt, extent = info1$extent, projection = info1$projection, 
              dimension = c(512, 512) * 2)
ximage::ximage(im, extent = info1$extent)
m <- do.call(cbind, maps::map(plot = F)[1:2])
lines(reproj::reproj_xy(m, info1$projstring), lwd = 5, col = "yellow")

image

{whatarelief} and {ximage} are two experiments for reading and visualizing, helpers on top of the GDAL and base R graphics to simplify work, they aren't finalized yet (gdalio was an earlier experiment)

Another example, get the last 21 days of images and see that we can get a quickview mosaic of them all on the same projection (we just get the extent of all of them, they are all the same projection so it's relatively easy).

library(raadfiles)
files <- polarview_files(type = "jpeg")
afile <- files$fullname[seq(nrow(files) - 20, nrow(files))]
library(whatarelief)
vrt <- character(length(afile))
exts <- matrix(0, nrow = length(vrt), ncol = 4)
for (i in seq_along(afile)) {
  tarball <- raadfiles:::polarview_jpeg_tarball(afile[i])

  info1 <- vapour::vapour_raster_info(glue::glue("/vsitar/{tarball}/{raadfiles:::polarview_tifname(tarball)}"))

  vrt0 <- vapour::vapour_vrt(afile[i], extent  = info1$extent, 
                          projection = info1$projection)

  print(info1$projstring)
  exts[i, ] <- info1$extent
   vrt[i] <- vrt0  
}

im <- imagery(source = vrt, extent = c(range(exts[,1:2]), range(exts[,3:4])), projection = info1$projection, 
              dimension = c(512, 512) * 2)
ximage::ximage(im, extent = c(range(exts[,1:2]), range(exts[,3:4])), asp = 1)
m <- do.call(cbind, maps::map(plot = F)[1:2])
lines(reproj::reproj_xy(m, info1$projstring), lwd = 1, col = "yellow")

image