search

lauraduncanson / icesat2_boreal

Biomass modeling and mapping of forest biomass in the boreal using NASA's ICESat-2
13 stars 6 forks source link

NAFlag isn't passed to WriteRaster #61

Open wildintellect opened 1 year ago

wildintellect commented 1 year ago

The Boreal Biomass tiles have Nodata set as nan but that should have been -9999, this makes the file metadata wrong, both for the Nodata, and the statistics. This was discovered by ORNL DAAC during submission.

Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  Description = lyr1
  Min=2.820 Max=184.518 
  Minimum=2.820, Maximum=184.518, Mean=-9999.000, StdDev=-9999.000
  NoData Value=nan
  Overviews: 1500x1500, 750x750, 375x375
  Metadata:
    STATISTICS_MAXIMUM=184.51815795898
    STATISTICS_MEAN=-9999
    STATISTICS_MINIMUM=2.8204426765442
    STATISTICS_STDDEV=-9999
Band 2 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = sd
  Min=0.414 Max=78.668 
  Minimum=0.414, Maximum=78.668, Mean=-9999.000, StdDev=-9999.000
  NoData Value=nan
  Overviews: 1500x1500, 750x750, 375x375
  Metadata:
    STATISTICS_MAXIMUM=78.668167114258
    STATISTICS_MEAN=-9999
    STATISTICS_MINIMUM=0.41386494040489
    STATISTICS_STDDEV=-9999

But this is an error where the NAflag was omitted from the writeRaster function in the code. https://github.com/lauraduncanson/icesat2_boreal/blob/b309d216af06d97be41469b1428b2fc984d16c84/lib/mapBoreal_simple.R#L1052-L1060

Should be:

writeRaster(out_map, filename=out_cog_fn, filetype="COG", NAflag=NAflag, gdal=c("COMPRESS=LZW", overwrite=TRUE, gdal=c("COMPRESS=LZW", "OVERVIEW_RESAMPLING=AVERAGE")))
wildintellect commented 1 year ago

To be clear we can also remove lines 1052 and 1054 that are reclassifying -9999 to NA.

wildintellect commented 1 year ago

Looking at this more in depth, based on the question from ORNL DAAC. Is nan or -9999.0 the desired NA value? Either way the stats in the file header are wrong (Mean and STDEV).

But we should agree on which value to use for NA before fixing the layer names and the stats (and figuring out what caused the bad stats to start, could be a bug in the WriteRaster).

pahbs commented 1 year ago 
# Here is what the code SHOULD look like before the writeRaster call
#

# remove some lines
#out_map <- subst(out_map, -9999, NA) 
out_sd <- app(out_map_all, sd) 
# out_sd <- subst(out_sd, -9999, NA) 
out_map <- c(out_map, out_sd) 

# To get our ndv recognized
NAflag(out_map) <- -9999

 tifoptions <- c("COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=6", "OVERVIEW_RESAMPLING=AVERAGE") 
 writeRaster(out_map, filename=out_cog_fn, filetype="COG", gdal=c("COMPRESS=LZW", overwrite=TRUE, gdal=c("COMPRESS=LZW", "OVERVIEW_RESAMPLING=AVERAGE")))