0

I am trying to view the following Drought Code file downloaded from the Zenodo server: [https://zenodo.org/record/3540926/files/drought_code_2017.nc?download=1], but the parts of Canada and the US are not visible in QGIS when I try to add the netCDF4 file as a raster layer (and mesh layer).

Screenshot: Mesh Layer import of the netCDF4 file

I tried the following things, but they did not work:

  • Changing the CRS to EPSG:4326 - WSG 84
  • Longitudes are in the range 0-360 in this netCDF4 file. I tried modifying it to set values in the range of -180 to 180
  • As recommended in the post Mapping a NetCDF in QGIS in the correct location on Earth, I tried running the command gdal_translate -a_srs EPSG:4326 -of netCDF drought_code_2017.nc drought_code_2017_4326.nc, but this did not work as well.

I am a student, new to QGIS.

[Added information-1] I tried retrieving values for near Forest McMurray in Alberta using Python and got the drought-code value:

>>> from netCDF4 import Dataset
>>> import numpy as np
>>> from datetime import datetime
>>> data = Dataset('drought_code_2017.nc')
>>> lat=56.752306
>>> lon=-111.476139 # Location in Alberta, near forest McMurray
>>> date = datetime.strptime("05/06/2017","%m/%d/%Y")
>>> day = (date.timetuple().tm_yday)-1 # It expects the day of year in number.
>>> # A sample subroutine to find the index of the nearest latitude and longitude from the netCDF File
>>> def find_nearest_value_index(array, value):
...     array = np.asarray(array)
...     idx = (np.abs(array - value)).argmin()
...     return idx
... 
>>> lon = 360 + lon if lon < 0 else lon # Normalizing the input longitude to the longitude range. The longitude in the dataset is represented in scale of 0-360. Input is -180 to 180
>>> lat_idx = find_nearest_value_index(data["Latitude"][:],lat)
>>> lon_idx = find_nearest_value_index(data["Longitude"][:], lon)
>>> data["DC"][day, lat_idx, lon_idx].data.item()
42.89582824707031
2
  • A lot of data north of 30°N are simply missing data I think. Its not only north-America and Canada, but also Europe and Asia. Looks like there are simply no measurements available in those spots for whatever reason (I'm not familiar with this particular dataset). What makes you think the data is not in the correct location? Sep 14, 2022 at 11:12
  • On their site, they mentioned that this is a Canadian Fire Weather Index dataset. Given a latitude and longitude in the region of Alberta (lat:56.752306, lon: -111.476139), I fetched the data using Python and was able to retrieve the data. I've included the Python commands in the post now. Sep 15, 2022 at 16:33

1 Answer 1

0

There simply is no data in that region on the map that you're showing and the value you found (42.89582824707031) is correct for that lat/lon/date.

I prefer using xarray, so I'm using that instead of netCDF4.

import datetime
import xarray as xr

# Open the dataset
ds = xr.open_dataset(r"/Users/hmcoerver/Downloads/drought_code_2017.nc", decode_coords="all")
# Create new longitude coordinates (to convert 0.360 -> -180.180)
new_longitude = np.where(ds.Longitude.values >= 180, ds.Longitude.values - 360, ds.Longitude.values)
# Assign the new coordinates
ds = ds.assign_coords({"Longitude": new_longitude}).sortby("Longitude")
# Select the variable.
da = ds["DC"]

# Define some point of interest.
lat = 56.752306
lon = -111.476139
doy = int(datetime.datetime.strptime("05/06/2017","%m/%d/%Y").strftime("%j"))

# Get the value at that point in time.
value_at_point = da.sel(Time = doy, Latitude = lat, Longitude = lon, method = "nearest").values

# Make a plot of the data.
fig = plt.figure(1)
ax = fig.gca()
da.sel(Time = doy).plot(ax = ax)
ax.scatter([lon], [lat], c = "tab:red")
ax.set_title(f"Day-of-year: {doy}, Value: {value_at_point:.4f}")
ax.set_facecolor("lightgray")
ax.grid()

You can see the value of 42.89582824707031 you found at the selected location (indicated by the red dot). Value of 42.8958 at the red dot.

The map you see in QGIS, is the first band of the dataset (i.e. at Day 1 of 2017), there is no valid data for (most of) Canada at that time:

doy = 1
fig = plt.figure(2)
ax = fig.gca()
da.sel(Time = doy).plot(ax = ax)
ax.scatter([lon], [lat], c = "tab:red")
ax.set_title(f"Day-of-year: {doy}, Value: {value_at_point:.2f}")
ax.set_facecolor("lightgray")
ax.grid()

There is no data at the selected point

We can also see that when plotting a time series at the selected location:

fig = plt.figure(3)
ax = fig.gca()
ts = da.sel(Latitude = lat, Longitude = lon, method="nearest").values
ax.plot(ts, label = "DC")
ax.plot(np.where(np.isnan(ts), 1, np.nan), label = "No Data")
ax.legend()
ax.set_title(f"lat: {lat:.3f}, lon: {lon:.3f}")
ax.set_facecolor("lightgray")
ax.grid()
ax.set_xlabel("Day-of-the-year")
ax.set_ylabel(f"{da.long_name}")

enter image description here

1
  • This was very detailed and helpful. Any suggestions on how to view the given day's data in QGIS? Sep 17, 2022 at 4:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.