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. For future reference, documentation on GDAL GRIB file unit conversion can be found here: <www.gdal.org/frmt_grib.html> With this information, I now realize that the previously mentioned 'GDAL values' were indeed my correct temperatures. However, I am still having trouble getting the correct latitude and longitude co-ordinates for each of these temperature values. I have trying to get the correct output for a single message/band for quite a while and 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]
All of the temperature data seems reasonable. However, the latitudes and longitudes don't seem to be correct. Every latitude-longitude pair seems to be tightly centered around 25 lat, -95 long. The problem with this is that the GRIB file I am using covers the entire continental U.S. and has a much wider range in lat-long values than the ones I am getting. 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