1

I have issues extracting geometrical shaped data from a GRIB raster of soil moisture. The file is for example https://gimms.gsfc.nasa.gov/SMOS/jbolten/FAS/L03/20100113_20100115.as1.grib

Using either GDAL or Rasterio I can read the file and display it here. enter image description here

Its CRS is CRS({'proj': 'stere', 'lat_0': 90, 'lat_ts': 60, 'lon_0': 1, 'k': 1, 'x_0': 0, 'y_0': 0, 'a': 6371200, 'b': 6371200, 'units': 'm', 'no_defs': True})

or with gdalinfo it reads

Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: 20141227_20141229.as1.grib
Size is 1440, 600
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371200,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",1],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (662126.628862303332426,59159218.394239947199821)
Pixel Size = (25000.000000000000000,-25000.000000000000000)
Corner Coordinates:
Upper Left  (  662126.629,59159218.394) (179d38'28.48"W, 67d16'31.51"S)
Lower Left  (  662126.629,44159218.394) (179d51'32.52"W, 59d52' 0.26"S)
Upper Right (36662126.629,59159218.394) (149d12'45.70"E, 70d36'45.13"S)
Lower Right (36662126.629,44159218.394) (141d17'58.79"E, 66d35'39.62"S)
Center      (18662126.629,51659218.394) (161d 8'15.06"E, 65d34'25.05"S)
Band 1 Block=1440x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC (Ground or water surface)
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=undefined [-]
    GRIB_ELEMENT=var255
    GRIB_FORECAST_SECONDS=86400 sec
    GRIB_REF_TIME=  1419638400 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[-]
    GRIB_VALID_TIME=  1419724800 sec UTC

The issue occurs if I want to extract a region from the raster, using a predefined shape. enter image description here

Using rasterio.warp.transform_geom with source CRS EPSG:4326 and dest CRS of the GRIB file above, I end up having a huge shape and of course the new one does not intersect with the raster. enter image description here

Could anyone show me how to do it right with this GRIB format? I have done such operation without any problem on netcdf, hdf, and geotiff files.

2

That data is definitely not in polar stereographic, that CRS has been misapplied - I had a stab in the dark, and assigned it longitude/latitude and bounded at 60S and 90N (using R):

library(raster)
r <- raster("~/20100113_20100115.as1.grib")
#plot(r) ## wrong
#plot(extent(r), add = T)

extent(r) <- extent(-180, 180, -60, 90)
projection(r) <- "+init=epsg:4326"
plot(r)
maps::map(add = T)

At the command line you can use gdal_translate arguments [-a_srs srs_def] [-a_ullr ulx uly lrx lry] to do the same, but - I would certainly want to know why this has occurred. You don't want to have to fix stuff like this ad hoc downstream.

Grib is a domain-specific format, and to have one created like this arbitrarily is weird to me - and it being done wrong is really going to cause problems for someone.

enter image description here

  • 1
    WYSIWYG turns out to be applied here. I am going to email the product admin if they ever want to support the users better. – Pateheo Oct 25 '18 at 9:50
  • Yeah, I was thinking there was a chance it might be GDAL for some reason loading the GRIB definition incorrectly, but I just opened it with NDFD DEGRIB and it's still showing that incorrect projection. Please contact the data admin and let them know about this issue! – Thomas Oct 25 '18 at 22:38
0

I found the correct answer very helpful, but I needed it in Python instead of R. So I've attached my code that (I believe) does the same thing as what @Thomas posted but in Python. The only Python library it uses is GDAL.

translate_options = gdal.TranslateOptions(format='GTiff', outputSRS='EPSG: 4326', outputBounds=[-179, 90, 180, -60])
translate_output = gdal.Translate(destName=output_filename, srcDS=src_filename, options=translate_options)

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.