AVHRR L1b GAC

[mdsumner]

Would like to arrange access to these products:

https://www.avl.class.noaa.gov/saa/products/search?datatype_family=AVHRR

The target is (see screenshot below)

• GAC
• Either descending or ascending
• any receiving station
• any time period
• intersection with the Antarctic coast, particularly Mawson

The files look like

• "NSS.GHRR.ND.D94182.S1520.E1650.B1625960.GC" (1994 example)
• "NSS.GHRR.TN.D79010.S1919.E2112.B0126061.GC" or "NSS.GHRR.TN.D79010.S1615.E1755.B0125859.WI" (1979 examples)

Don't yet know how to parse out target files for a given region - but maybe there's not so much data that it matters - we have 98 files that cross near Mawson in 1979 and they total 60Mb (AFraser found by search with the website).

I don't know what GC and WI stand for or any of the other file tokens except "D94182" is "1994, day 182" .

This is the GDAL driver: https://gdal.org/frmt_l1b.html the files work well with gdalwarp, by specifying a target region and projection with the "-geoloc" option.

Application is fast ice detection (manual) from visual inspect of band 4 or 5.

gdalinfo -nogcp yields:

Driver: L1B/NOAA Polar Orbiter Level 1b Data Set
Files: Mawson1994/GCdata/NSS.GHRR.ND.D94182.S1520.E1650.B1625960.GC
Size is 409, 10764
Coordinate System is `'
Metadata:
  DATASET_NAME=NSS.GHRR.ND.D94182.S1520.E1650.B1625960.GC
  DATA_TYPE=AVHRR GAC
  LOCATION=Descending
  PROCESSING_CENTER=NOAA/NESDIS - Suitland, Maryland, USA
  REVOLUTION=16259
  SATELLITE=NOAA-12(D)
  SOURCE=Fairbanks, Alaska, USA (formerly Gilmore Creek)
  START=year: 1994, day: 182, millisecond: 55247146
  STOP=year: 2000, day: 0, millisecond: 0
Subdatasets:
  SUBDATASET_1_NAME=L1B_SOLAR_ZENITH_ANGLES:"Mawson1994/GCdata/NSS.GHRR.ND.D94182.S1520.E1650.B1625960.GC"
  SUBDATASET_1_DESC=Solar zenith angles
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  SRS=GEOGCS["WGS 72",DATUM["WGS_1972",SPHEROID["WGS 72",6378135,298.26,AUTHORITY["EPSG",7043]],TOWGS84[0,0,4.5,0,0,0.554,0.2263],AUTHORITY["EPSG",6322]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",9108]],AUTHORITY["EPSG",4322]]
  X_BAND=1
  X_DATASET=L1BGCPS_INTERPOL:"Mawson1994/GCdata/NSS.GHRR.ND.D94182.S1520.E1650.B1625960.GC"
  Y_BAND=2
  Y_DATASET=L1BGCPS_INTERPOL:"Mawson1994/GCdata/NSS.GHRR.ND.D94182.S1520.E1650.B1625960.GC"
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,10764.0)
Upper Right (  409.0,    0.0)
Lower Right (  409.0,10764.0)
Center      (  204.5, 5382.0)
Band 1 Block=409x1 Type=UInt16, ColorInterp=Undefined
  Description = AVHRR Channel 1:  0.58  micrometers -- 0.68 micrometers
Band 2 Block=409x1 Type=UInt16, ColorInterp=Undefined
  Description = AVHRR Channel 2:  0.725 micrometers -- 1.10 micrometers
Band 3 Block=409x1 Type=UInt16, ColorInterp=Undefined
  Description = AVHRR Channel 3:  3.55  micrometers -- 3.93 micrometers
Band 4 Block=409x1 Type=UInt16, ColorInterp=Undefined
  Description = AVHRR Channel 4:  10.3  micrometers -- 11.3 micrometers
Band 5 Block=409x1 Type=UInt16, ColorInterp=Undefined
  Description = AVHRR Channel 5:  11.5  micrometers -- 12.5 micrometers

[raymondben]

That search interface is unfortunately a bit too intractable for our current methods.

[mdsumner]

Files

I put the current collection into

data_local/www.avl.class.noaa.gov/saa/products/AVHRR_L1b/

The AVHRR_L1b is a dir name made up by me, it contains 1978/ 1979/ and 1994/ based on .D%y, i.e. NSS.GHRR.ND.D94182 is 1994 jday 182.

The extensions are WI or GC, which are the receiving stations (you can read this stuff in the gdalinfo output).

Getting more

afaik we must use the website search to get a file listing, and this can only be done one year at a time - I think it's about 32Gb per year for all so a fun rainy day task.

reading, warping

This code is enough to build the gdalwarp command for a target raster. I believe GDAL 2.3.2 is better than 2.2.3 for this task - and possibly it's an unfinished facility for some problems. The GCPs are a subsample of the grid itself.

GDAL uses the internal GCPs (also obtainable with vapour::vapour_raster_gcp in pixel, line, x, y, z form.

#' Generate gdalwarp command 
#'
#' A raster object 'x' is used to input the extent, projection (CRS), and
#' resolution (in metres) of the target grid. 
#' 
#' The output is a templated string, which can be sprintf'ed to apply
#' input and output file name. 
#' 
#' Resolution is currently hardcoded at 1000x1000m, so that the input can
#' be any general grid that we are familiar with (WIP). 
#' 
#' @param x raster object, the reprojection target
#' @param resampling method, default is bilinear
#'
#' @return
#' @export
#'
#' @examples
generate_gdalwarp_l1b <- function(x, resampling = "bilinear") {
  p4 <- projection(x)
  te <- c(xmin(x), ymin(x), xmax(x), ymax(x))
  tr <- res(x)
  ## input output resampling t_srs te tr
  paste("gdalwarp %s %s", sprintf("-geoloc -r %s -t_srs '%s' -te %f %f %f %f -tr %f %f -co COMPRESS=LZW -dstnodata 0",
                                  resampling,
                                  p4, 
                                  te[1], te[2], te[3], te[4],
                                  tr[1], tr[2]
  ))
}
#' Run GDAL warp 
#'
#' All bands in the input are written to the output. The file is 
#' compressed using LZW. Warping is done via geolocation arrays, which 
#' GDAL just handles on its own. 
#' @param input NOAA AVHRR filename
#' @param target raster object (composed of extent, crs, and resolution)
#' @param output output file name (*.tif) which is GeoTIFF
#' @param resampling resampling method, default is 'bilinear' (gdalwarp default is 'near')
#'
#' @return raster object from 'output' file
#' @export
#'
gdalwarp <- function(input, target, output, resampling = "bilinear", overwrite = FALSE) {
  cmd <- sprintf(generate_gdalwarp_l1b(target, resampling = resampling), input, output)
  if (overwrite) {
    cmd <- paste(cmd, "-overwrite")
  } else {
    if (file.exists(output)) stop("output file exists, use 'overwrite = TRUE' or delete it")
  }
  tst <- system(cmd)
  brick(output)
}