OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.95k stars 2.57k forks source link

Area and Temp bands show no data in GOES FDCF fire data #1811

Closed schwehr closed 5 years ago

schwehr commented 5 years ago

This may be already fixed in head. I am currently investigating. I've been using an older version from debian testing. There could be an error in my usage of gdal, in gdal reading the netcdf, or an error in the convention inside of the GOES files. I checked through ~20K files with GDAL and got no valid pixels in Area or Temp.

gdalinfo --version
GDAL 2.3.2, released 2018/09/21

I am able to see data in the Mask, Power, and DQF bands. I get a

gsutil cp gs://gcp-public-data-goes-16/ABI-L2-FDCF/2019/200/17/OR_ABI-L2-FDCF-M6_G16_s20192001740408_e20192001750116_c20192001750220.nc .

gdalinfo -stats -mm -hist -checksum NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Mask | egrep 'bucket|Checksum|NoData|Fill|Type|Mean'
  Mask#_FillValue=-99
Band 1 Block=5424x1 Type=Int16, ColorInterp=Undefined
  Minimum=10.000, Maximum=245.000, Mean=136.997, StdDev=44.372
  256 buckets from 9.53922 to 245.461:
  Checksum=22221
  NoData Value=-99
  Unit Type: 1
    _FillValue=-99

gdalinfo -stats -mm -hist -checksum NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Area | egrep 'bucket|Checksum|NoData|Fill|Type|Mean'
  Area#_FillValue=-1
ERROR 1: NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Area, band 1: Failed to compute min/max, no valid pixels found in sampling.
ERROR 1: NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Area, band 1: Failed to compute statistics, no valid pixels found in sampling.
ERROR 1: NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Area, band 1: Failed to compute statistics, no valid pixels found in sampling.
Band 1 Block=5424x1 Type=Int16, ColorInterp=Undefined
  Checksum=5888
  NoData Value=-1
  Unit Type: m2
    _FillValue=-1

Checking via python, it looks like there are no pixels other than -1 coming back from the file.

ipython3
from osgeo import gdal
src = gdal.Open('NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Area')
src.RasterCount
Out[58]: 1
band = src.GetRasterBand(1)
data = band.ReadAsArray()
data
Out[61]: 
array([[-1, -1, -1, ..., -1, -1, -1],
       [-1, -1, -1, ..., -1, -1, -1],
       [-1, -1, -1, ..., -1, -1, -1],
       ...,
       [-1, -1, -1, ..., -1, -1, -1],
       [-1, -1, -1, ..., -1, -1, -1],
       [-1, -1, -1, ..., -1, -1, -1]], dtype=int16)

import numpy as np
np.where( data > -1)
Out[62]: (array([], dtype=int64), array([], dtype=int64))

Double checking with xarray, I do see data:

import xarray
nc = xarray.open_dataset('OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc')
nc
nc.Area
nc.Area.data
data = nc.Area.data

import numpy as np
np.count_nonzero(~np.isnan(data))

places = list(map(tuple, np.where(~np.isnan(data))))
len(places[0]) # 222

for i in range (len(places[0])):
    x,y = places[0][i], places[1][i]
    print (x, y, data[x, y])

[SNIP]
3308 5234 4000.0
3308 5235 4000.0
3309 5234 4000.0
3309 5235 4000.0
3310 5234 4000.0
3310 5235 4000.0
3311 5234 4000.0
3750 409 546478.06
3750 410 180110.23

A sample visualization in Earth Engine so folks have a sense of what it looks like:

goes-fdcf-demo
schwehr commented 5 years ago

Setting HONOUR_VALID_RANGE to false appears to make it work.

A bit more info. The gdalinfo that I'm trying with reports:

gdalinfo --format netcdf
Format Details:
  Short Name: netCDF
  Long Name: Network Common Data Format
  Supports: Raster
  Supports: Vector
  Extension: nc
  Help Topic: frmt_netcdf.html
  Supports: Subdatasets
  Supports: Open() - Open existing dataset.
  Supports: Create() - Create writable dataset.
  Supports: CreateCopy() - Create dataset by copying another.
  Creation Field Datatypes: Integer Integer64 Real String Date DateTime

<CreationOptionList>
[SNIP]
</CreationOptionList>

<LayerCreationOptionList>
[SNIP]
</LayerCreationOptionList>

<OpenOptionList>
  <Option name="HONOUR_VALID_RANGE" type="boolean" description="Whether to set to nodata pixel values outside of the validity range" default="YES" />
</OpenOptionList>

  Other metadata items:
    GDAL_HAS_HDF4=YES
    GDAL_HAS_HDF5=YES
    NETCDF_CONVENTIONS=CF-1.5
    NETCDF_HAS_NC2=YES
    NETCDF_HAS_NC4=YES
    NETCDF_VERSION=4.6.2 of Nov 20 2018 06:04:35 $

gdalinfo reports this about the file:

gdalinfo NETCDF:OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc:Area | egrep -i 'convention|1\.7|vocab' | grep NC
  NC_GLOBAL#Conventions=CF-1.7
  NC_GLOBAL#keywords_vocabulary=NASA Global Change Master Directory (GCMD) Earth Science Keywords, Version 7.0.0.0.0
  NC_GLOBAL#Metadata_Conventions=Unidata Dataset Discovery v1.0
  NC_GLOBAL#standard_name_vocabulary=CF Standard Name Table (v35, 20 July 2016)
schwehr commented 5 years ago

The files have invalid ranges for Area and Temp:

ncdump OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc | grep valid_range *.ncdump.txt
                Area:valid_range = 0s, -6s ;  <-- OUCH
                Temp:valid_range = 0s, -6s ;  <-- OUCH
                Mask:valid_range = 0s, 245s ;
                Power:valid_range = 75.f, 50000.f ;
                DQF:valid_range = 0b, 5b ;
                minimum_fire_temperature:valid_range = 600.f, 1200.f ;
                maximum_fire_temperature:valid_range = 600.f, 1200.f ;
                mean_fire_temperature:valid_range = 600.f, 1200.f ;
                minimum_fire_area:valid_range = 4000.f, 4000000.f ;
                maximum_fire_area:valid_range = 4000.f, 4000000.f ;
                mean_fire_area:valid_range = 4000.f, 4000000.f ;
                minimum_fire_radiative_power:valid_range = 75.f, 50000.f ;
                maximum_fire_radiative_power:valid_range = 75.f, 50000.f ;
                mean_fire_radiative_power:valid_range = 75.f, 50000.f ;
                percent_uncorrectable_GRB_errors:valid_range = 0.f, 1.f ;
                percent_uncorrectable_L0_errors:valid_range = 0.f, 1.f ;
schwehr commented 5 years ago

Followup:

A python demonstration of HONOUR_VALID_RANGE working:

  def testHonourValidRange(self):
    # Area has an invalid valid_range of 0:-6
    filename = 'OR_ABI-L2-FDCF-M6_G17_s20192011300341_e20192011309408_c20192011309501.nc'
    filepath = 'NETCDF:' + self.getTestFilePath(filename) + ':Area'

    # ERROR 1: band 1:
    # Failed to compute statistics, no valid pixels found in sampling.
    with gcore_util.ErrorHandler('CPLQuietErrorHandler'):
      src = gdal.OpenEx(filepath)
      band = src.GetRasterBand(1)
      stats = band.GetStatistics(False, True)
      self.assertEqual([0.0, 0.0, 0.0, 0.0], stats)
      src = None

    src = gdal.OpenEx(filepath, open_options=['HONOUR_VALID_RANGE=FALSE'])
    band = src.GetRasterBand(1)
    stats = band.GetStatistics(False, True)
    self.assertEqual([0.0, 8896.0, 53.08108108108107, 625.4866665265016], stats)
    src = None

    src = gdal.OpenEx(filepath, open_options=['HONOUR_VALID_RANGE=TRUE'])
    band = src.GetRasterBand(1)
    stats = band.GetStatistics(False, True)
    self.assertEqual([0.0, 0.0, 0.0, 0.0], stats)
    src = None

The issue is known to the processing team. Their readme's are in the process of being updated, so don't have all the details yet

https://www.ncdc.noaa.gov/data-access/satellite-data/goes-r-series-satellites#FDC

schwehr commented 5 years ago

It's possible to fix the netcdf GOES files like this (and for Temp and Power):

    dataset = netCDF4.Dataset(dst, 'r+')
    area = dataset['Area']
    area.setncattr('valid_range', numpy.array([-32768, 32767], dtype=numpy.int16))
    dataset.close()
    del dataset

Or remove the valid_range.