TransWikia.com

GDAL: wrong corner coordinates from KNMI grib file

Geographic Information Systems Asked by Carpenter on April 17, 2021

I’m trying to parse KNMI grib files to extract weather forecasts for The Netherlands. I use gdal.netcore which is a .NET wrapper for gdal 3.2.0. On the website of KNMI they say the following:

  • Rekenrooster: Lambertprojectie
  • Product resolutie: 0,037° west-oost en 0,023° noord-zuid (300×300 punten)
  • Grid representatie: Regulier Lat/Lon
  • Gebiedsdefinitie: noord 55,877° N.B. zuid 49,0° N.B. / oost 11,063° O.L. en west 0,0° W.L.

However when I run

var info = Gdal.GDALInfo(dataset, new GDALInfoOptions(null));
Console.WriteLine(info);

I get the following information:

Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: /home/HA40_N25_202103190600_00000_GB
Size is 300, 300
Coordinate System is:
GEOGCRS["Coordinate System imported from GRIB file",
    DATUM["unnamed",
        ELLIPSOID["Sphere",6367470,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (-0.018500000000000,55.888500000000001)
Pixel Size = (0.037000000000000,-0.023000000000000)
Corner Coordinates:
Upper Left  (  -0.0185000,  55.8885000) (  0d 1' 6.60"W, 55d53'18.60"N)
Lower Left  (  -0.0185000,  48.9885000) (  0d 1' 6.60"W, 48d59'18.60"N)
Upper Right (  11.0815000,  55.8885000) ( 11d 4'53.40"E, 55d53'18.60"N)
Lower Right (  11.0815000,  48.9885000) ( 11d 4'53.40"E, 48d59'18.60"N)
Center      (   5.5315000,  52.4385000) (  5d31'53.40"E, 52d26'18.60"N)

This is quite close (55.8885 vs 55.877) but nonetheless inaccurate. Even if I try to project it, it still shows the same coordinates.

GdalBase.ConfigureAll();
using var dataset = Gdal.Open(@"D:HA40_N25_202103190600_00000_GB", Access.GA_ReadOnly);
var info = Gdal.GDALInfo(dataset, new GDALInfoOptions(null));

Console.WriteLine(info);

double[] adfGeoTransform = new double[6];
dataset.GetGeoTransform(adfGeoTransform);

var upperX = adfGeoTransform[0];
var upperY = adfGeoTransform[3];

var lowerX = upperX + (dataset.RasterXSize * adfGeoTransform[1]);
var lowerY = upperY + (dataset.RasterYSize * adfGeoTransform[5]);

var projection = dataset.GetProjection();
var source = new OSGeo.OSR.SpatialReference(null);
source.ImportFromWkt(ref projection);

var target = new OSGeo.OSR.SpatialReference(null);
target.ImportFromEPSG(4326);

var transform = new OSGeo.OSR.CoordinateTransformation(source, target);
double[] data = new double[3];
transform.TransformPoint(data, upperX, upperY, 0);

Console.WriteLine(data[0].ToString() + " " + data[1].ToString());

So what am I doing wrong here?

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP