4

I have a netCDF file with spatial data ranging from x:528500 573500 y:350500 399500 on EPSG: 31287. If I open and explore the file with Panoply and show the data on a map all works like it should.

But if I open it in QGIS I got a strange left shift of 360m in the x coordinate.

So I looked for a potential issue in the netCDF metadata with gdalinfo, and found, that the spatial layer extent causes the problem. I have a python program that automatically produces netCDF files using the NETCDF4-python library (Projection and other attributes are assigned regarding to the NetCDF Coordinate Convention Proposals).

Does anybody know why the extent is not exactly the same as the netCDF files x/y coordinates?

Here goes my gdalinfo of the netCDF's data variable:

Driver: netCDF/Network Common Data Format
Files: testdata.nc
Size is 46, 50
Coordinate System is:
PROJCS["MGI / Austria Lambert",
    GEOGCS["MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4312"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",49],
    PARAMETER["standard_parallel_2",46],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",400000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","31287"]]
Origin = (527640.000000000000000,400000.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
  Lambert_Conformal_Conic#false_easting=400000
  Lambert_Conformal_Conic#false_northing=400000
  Lambert_Conformal_Conic#GeoTransform=111000 69.64118155840001 0 572000 0 -89.4972135177 
  Lambert_Conformal_Conic#grid_mapping_name=lambert_conformal_conic
  Lambert_Conformal_Conic#inverse_flattening=299.1528128
  Lambert_Conformal_Conic#latitude_of_projection_origin=47.5
  Lambert_Conformal_Conic#longitude_of_central_meridian=13.33333333333333
  Lambert_Conformal_Conic#longitude_of_prime_meridian=0
  Lambert_Conformal_Conic#long_name=CRS definition
  Lambert_Conformal_Conic#semi_major_axis=6377387.155
  Lambert_Conformal_Conic#spatial_ref=PROJCS["MGI / Austria Lambert",GEOGCS["MGI",DATUM["Militar_Geographische_Institute",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],AUTHORITY["EPSG","6312"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4312"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",46],PARAMETER["latitude_of_origin",47.5],PARAMETER["central_meridian",13.33333333333333],PARAMETER["false_easting",400000],PARAMETER["false_northing",400000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","31287"]]
  Lambert_Conformal_Conic#standard_parallel={49,46}
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={1,6}
  NETCDF_DIM_time_VALUES=7706.5
  tas#cell_methods=time: mean
  tas#coordinates=x y
  tas#grid_mapping=Lambert_Conformal_Conic
  tas#long_name=Daily Mean Near-Surface Air Temperature
  tas#standard_name=air_temperature
  tas#units=Celsius
  tas#_FillValue=1e+020
  time#units=days since 1949-12-01 00:00:00
  x#axis=X
  x#long_name=x coordinate of projection
  x#standard_name=projection_x_coordinate
  x#units=meters
  y#axis=Y
  y#long_name=y coordinate of projection
  y#standard_name=projection_y_coordinate
  y#units=meters
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
  X_BAND=1
  X_DATASET=NETCDF:"testdata.nc":x
  Y_BAND=1
  Y_DATASET=NETCDF:"testdata.nc":y
Corner Coordinates:
Upper Left  (  527640.000,  400000.000) ( 15d 1'40.64"E, 47d29'14.91"N)
Lower Left  (  527640.000,  350000.000) ( 15d 0'48.97"E, 47d 2'15.52"N)
Upper Right (  573640.000,  400000.000) ( 15d38'18.12"E, 47d28'36.56"N)
Lower Right (  573640.000,  350000.000) ( 15d37' 7.86"E, 47d 1'37.50"N)
Center      (  550640.000,  375000.000) ( 15d19'28.88"E, 47d15'27.59"N)
Band 1 Block=46x1 Type=Float32, ColorInterp=Undefined
  NoData Value=1.0000000200408773e+020
  Unit Type: Celsius
  Metadata:
    cell_methods=time: mean
    coordinates=x y
    grid_mapping=Lambert_Conformal_Conic
    long_name=Daily Mean Near-Surface Air Temperature
    NETCDF_DIM_time=7706.5
    NETCDF_VARNAME=tas
    standard_name=air_temperature
    units=Celsius
    _FillValue=1e+020

Regarding to the the pixel size of 1000x1000 and a dataset dimension of 46,50 the appropriate extent should be x:528000 574000 y:350000 400000, but is calculated/set to x:527640 573640 y:350000 400000 on file creation automatically.

A possible solution to solve the issue is to reprocess the data with gdalwarp and to set manually the extent with -te: gdalwarp -te 528000 350000 574000 400000 NETCDF:"testfile.nc:tas" 1_proj_final.nc -of netCDF.

But that's not the way I would like to do it, so I am very interested in, to find out, why the extend is not automatically set correctly on file creation! Does anybody have an idea which part causes the issue ?

HERE goes the testfile netCDF.

  • Are you sure this isn't an issue with WGS84 versus MGI? MGI is a few hundred meters off from WGS84. – mkennedy Jan 2 '17 at 0:21
  • @mkennedy this would not be exactly 360 meters, and only for the X coordinate. – AndreJ Jan 2 '17 at 7:39
  • 1
    I can reproduce the bug, and I think it is worth reporting at trac.osgeo.org/gdal/report and the gdal-dev mailing list. The file contains 1-D geolocation arrays, which might not be used correctly by GDAL. Maybe it is trying to subtract 360 degrees instead of using meters as units. – AndreJ Jan 2 '17 at 8:16
  • Thank you for helping. Using the netCDF with a spatial reference, the extent have to be in 1-D arrays. I tried it with 2-D arrays but that doestn work, neither in QGIS nor in Panoply. – Manuel Jan 5 '17 at 15:04
  • @mkennedy Hmm maybe the SRS, which is in a WGS84 Datum causes the problem, but i thought if I create the file with a MGI Datum all other features depends on this. Is it posssible to change the SRS in the netCDF file on file creation ? – Manuel Jan 5 '17 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.