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

polarview JPEG images #126

Open
mdsumner opened this issue Jan 20, 2021 · 1 comment
Open

polarview JPEG images #126

mdsumner opened this issue Jan 20, 2021 · 1 comment

Comments

@mdsumner
Copy link
Member

mdsumner commented Jan 20, 2021

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
Copy link
Member Author

mdsumner commented Aug 30, 2022

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

@mdsumner mdsumner transferred this issue from AustralianAntarcticDivision/raadfiles Aug 30, 2022
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