0

I want to determine the lat/long coordinates of my raster layer. I calculated them using the geotransform of the dataset.

When I run the script on my file, it produces longitudes that are out of range:

(0.0, 0.0,  '->', 179.88, 90.12)
(0, 751,    '->', 179.88, -90.11999999999998) 
(1500, 751, '->', 539.88, -90.11999999999998) 
(1500, 0,   '->', 539.88, 90.12)

The longitude ranges from 180 to 540 (it's off by 360)

Here's the code I used:

import gdal, osr

gribdata = gdal.Open('CMC_glb_ABSV_ISBL_250_latlon.24x.24_2018041800_P000.grib2')
gribdata_gtrn = gribdata.GetGeoTransform()

gribdata_bbox_cells = (
    (0., 0.),
    (0, gribdata.RasterYSize),
    (gribdata.RasterXSize, gribdata.RasterYSize),
    (gribdata.RasterXSize, 0),
  )

geo_pts = []
for x, y in gribdata_bbox_cells:
    x2 = gribdata_gtrn[0] + gribdata_gtrn[1] * x + gribdata_gtrn[2] * y
    y2 = gribdata_gtrn[3] + gribdata_gtrn[4] * x + gribdata_gtrn[5] * y
    print (x, y, '->', x2, y2)

The file I'm using is this one, taken from the CMC global forecast.

The latitudes are correct (90 to -90), but howcome the longitude is 180 to 540?

stanri
  • 265
  • 3
  • 10
  • Your grib2 file is certainly not projected (you can check that with gribdata_srs.IsProjected()), so of course the coordinate transformation won't work – obchardon Apr 18 '18 at 14:15
  • So indeed, the raster is not projected. You will have to get the information about the projection system of this data. – obchardon Apr 18 '18 at 14:27
  • @obchardon you are correct, but the question still stands because even without the transform the longitudes are off by 360. I have edited and included an example without the transform. – stanri Apr 18 '18 at 14:39
  • 1
    But how can you expect the output coordinate to be right if you do not set the projection system, there is no error in this code ! – obchardon Apr 18 '18 at 15:18

0 Answers0