0

Suppose in Python, I have read in a netCDF4 file, as follows:

import netCDF4

ds = netCDF4.Dataset('fire_weather_index_2018.nc')

The aspects of the netCDF4 is as follows:

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    Conventions: CF-1.4
    created_by: R, packages ncdf4 and raster (version 3.0-7)
    date: 2019-10-16 00:07:39
    dimensions(sizes): Longitude(1440), Latitude(721), Time(365)
    variables(dimensions): int32 crs(), float64 Longitude(Longitude), float64 Latitude(Latitude), float64 Time(Time), float32 FWI(Time, Latitude, Longitude)
    groups:

When trying to access a particular datapoint in the file (for example, for day 200, at the coordinates of (27, 47), though this error occurs for every datapoint), as follows:

ds['FWI'][200, 27, 47]

It says the data is masked, but when I query directly ds['FWI'], there are plenty of datapoints. In fact, I can map it accordingly on a world map and all the data would show up.

How do I access / extract a particular data point for the 'FWI' variable in this case?

NicTam
  • 1
  • 1
    When you say coordinate, do you mean latitude/longitude? – snowman2 Nov 12 '21 at 14:36
  • @snowman2 yes, that's right – NicTam Nov 12 '21 at 14:55
  • This may help: https://gis.stackexchange.com/questions/358036/extracting-data-from-a-raster/358058#358058 – snowman2 Nov 12 '21 at 15:46
  • You might look in your Longitude(1440), Latitude(721) variables to see which index matches your coordinates. You may be getting a masked value at whatever coordinate happens to be stored at Latitude(27) Longitude(47) in those arrays. – Dave X Nov 15 '21 at 19:09
  • Based on the 721 x 1440, the data array is probably 0.25° resolution and you are indexing 27,47 in it. ds['FWI'][200, 27, 47] could well be at 27/72090-90=-86.625° 47/1440180-180=-174.125°, signs depending on the ordering of the Latitude, longitude variables. – Dave X Nov 16 '21 at 22:25
  • Is your dataset like those available from https://zenodo.org/record/3539654#.YajbsvHMLMI ? – Dave X Dec 02 '21 at 15:39

1 Answers1

0

You are getting a masked value because you are ignoring a level of abstraction and treating your coordinates as indices. The dataset probably holds a masked result for [27,47] at 83.25°N, 11.75°E, which is well north of Svalbard.

You need to use the indices that match your coordinates:

import netCDF4
import numpy as np

ds = netCDF4.Dataset('fire_weather_index_2018.nc')


print(ds['Latitude'][27])
print(ds['Longitude'][47])


jj = np.argwhere(ds['Latitude'][:]==27)[0][0]  # first matching lat
ii = np.argwhere(ds['Longitude'][:]==47)[0][0] # first matching lon


[ii,jj]


ds['FWI'][200, jj, ii]

If the coordinates don't have an exact match (e.g. 26.3°N,47.3°E), you can look for the nearest ones:

jj = np.argmin((ds['Latitude'][:]-27.3)**2)
ii = np.argmin((ds['Longitude'][:]-47.3)**2)

print([ii,jj])
print(ds['Latitude'][jj])
print(ds['Longitude'][ii])
print(ds['FWI'][200, jj, ii])

(I found some comparable data at https://zenodo.org/record/3539654#.YajbsvHMLMI )

Dave X
  • 1,639
  • 13
  • 25