I am attempting to use GDAL to read corresponding latitude, longitude, and temperature values from a GRIB file. I would like to store these data points in a 2-D list. I was able to do this using the NDFD GRIB Decoder. However, this tool required me to first create a CSV of the data. If possible, I would like to avoid creating any files because I will be reading a large number of GRIB files. Thus, I switched to GDAL.
Currently, I am able to use GDAL to iterate through each message/band of the GRIB file and view the metadata of each message/band. However, I don't understand how to extract the data from the message/band. I tried using ReadAsArray() and got a large 2-D list of mostly 9999's (No Data Value). I iterated through this list to do some quick stats on all the non-'No Data Values'.
dataset = gdal.Open('grib_file.grb', gdal.GA_ReadOnly)
message = dataset.GetRasterBand(1)
data_array = message.ReadAsArray()
num_list = []
for row in data_array:
for value in row:
if value < 9999:
num_list.append(value)
print("Count: " + str(len(num_list)))
print("Max: " + str(np.max(num_list)))
print("Min: " + str(np.min(num_list)))
print("Mean: " + str(statistics.mean(num_list)))
print("Standard Deviation: " + str(statistics.stdev(num_list)))
Comparing the resulting stats to a 'degribbed' CSV of the same message/band, I found that the number of non-9999 GDAL values in my 2-D list is equal the number of data rows in my CSV. So I assume each of these GDAL values correlates to 1 latitude-longitude-temperature data point. However, when looking at the max, min, and mean of the GDAL values, they don't match my CSV's values for latitude, longitude, or temperature. What do these GDAL values actually represent? And what is the best way for me to extract the latitude-longitude-temperature data points?
I am very new to working with GIS data, so I may have a fundamental misunderstanding as to how GDAL or GRIB files work.
Update
After some more research, I found that GDAL automatically converts all temperature values in GRIB files to Celsius by default. Additionally, the NDFD Decoder (aka 'degribber') converts to Fahrenheit and rounds to the nearest integer. My temperature values were correct, just using a different unit than I expected. For future reference, documentation on GDAL GRIB file unit conversion can be found here: https://www.gdal.org/frmt_grib.html.
However, I am still having trouble getting the correct latitude and longitude co-ordinates for each of these temperature values. Here is what I have right now:
import gdal
import numpy as np
import statistics
import osr
import math
# Open file
dataset = gdal.Open('E:/Downloads/YEUZ98_KWBN_201001011259.grb2', gdal.GA_ReadOnly)
message_count = dataset.RasterCount
x_size = dataset.RasterXSize
y_size = dataset.RasterYSize
# Preparing transformation
src_srs = osr.SpatialReference()
src_srs.ImportFromWkt(dataset.GetProjection())
tgt_srs = osr.SpatialReference()
tgt_srs.ImportFromEPSG(4326)
transform = osr.CoordinateTransformation(src_srs, tgt_srs)
# Parsing for valid data points
message = dataset.GetRasterBand(1)
data_array = message.ReadAsArray()
data_points = []
for row in range(y_size):
for col in range(x_size):
temperature = data_array[row][col]
if temperature != message.GetNoDataValue():
lat_long_point = transform.TransformPoint(row, col)
lat = lat_long_point[1]
long = lat_long_point[0]
data_points.append([lat, long, temperature])
# Display statistics for temperature
temperatures = [data_point[2] for data_point in data_points]
print("Count: " + str(len(temperatures)))
print("Max: " + str(np.max(temperatures)))
print("Min: " + str(np.min(temperatures)))
print("Mean: " + str(statistics.mean(temperatures)))
print("Standard Deviation: " + str(statistics.stdev(temperatures)))
# Show 1/20000 of the data points. Each data point holds a temperature and its corresponding lat/long
print("\nData Points:")
for i in range(math.floor(len(data_points) / 20000)):
print(data_points[i * 20000])
I get the following output:
Count: 368246
Max: 24.950006103515648
Min: -31.649999999999977
Mean: -4.05918937533062
Standard Deviation: 10.215615846529928
Data Points [Latitude, Longitude, Temperature]:
[25.000890299683032, -94.99952371153155, 5.550012207031273]
[25.00491913062724, -94.99888862379909, -25.549993896484352]
[25.001070152573444, -94.99860090057608, -0.04999389648435226]
[25.00069244683015, -94.99835283836573, 7.249993896484398]
[25.005575607284577, -94.99813446569522, -14.449987792968727]
[25.001942459037867, -94.99790629734547, -2.2500061035156023]
[25.00026077721243, -94.99768802823817, 12.249993896484398]
[25.00249102141579, -94.99747960735287, -6.649999999999977]
[25.00358815549422, -94.99726128117258, -9.449987792968727]
[25.008084618188125, -94.99704286927876, 3.8500000000000227]
[25.000404647896634, -94.99680491081958, 11.150018310546898]
[25.003300367512978, -94.99657660948411, -10.049993896484352]
[25.00069241108848, -94.99633853829565, 12.249993896484398]
[25.00553959660932, -94.99611016365309, -3.3499816894531023]
[25.00364208093872, -94.99586215125959, -3.8499816894531023]
[25.006762621284626, -94.995594121097, 6.150018310546898]
[25.00740110988653, -94.99527655572876, 13.350000000000023]
[25.003453198251332, -94.9948599463964, -0.04999389648435226]
Every latitude-longitude pair seems to be tightly centered around 25 lat, -95 long. However, the GRIB file I am using covers the entire continental U.S. What am I doing wrong?
Downloads for the GRIB file I am working with and the CSV for the 1st message/band of that GRIB can be found here: https://www.dropbox.com/sh/oiaq91jq27isbp8/AADlskNq68sC_dhb7J8GXfBaa?dl=0
Best Answer
I downloaded a free grib file from this link:
gribs2.gmn-usa.com/cgi-bin/weather_fetch.pl?parameter=wind&days=7®ion=Pacific&dataset=nww3
and it was opened with QGIS (it uses gdal and grib raster format is clearly supported); as it can be observed at following image:
In this particular case there are 54 bands and their statistics are different for each one of them.
Adapting your code for calculating statistics in each particular case (it was corroborated that 999 effectively corresponds to nodata values with band.GetNoDataValue()):
I got following results after running above code at Python Console of QGIS: