[GRASS-user] Import and reproject COSMO REA6 data set

[Helmut Kudrnovsky]

Pietro Zambelli wrote
> Dear all,
> 
> I'm not able to properly import the Cosmo regional reanalysis (COSMO
> REA6) data set.
> 
> http://reanalysis.meteo.uni-bonn.de/?COSMO-REA6
> 
> The original file is available in a GRIB format here:
> 
> ftp://ftp-cdc.dwd.de/pub/REA/COSMO_REA6/hourly/2D/T_2M/T_2M.2D.201501.grb.bz2
> 
> If I import the cosmo rea6 data set and display also the National
> statistical units for the EU, The two data sets don't match, see the
> attached figure.
> 
> What I'm doing wrong?
> How should I properly importing the COSMO REA6 data set?
> Any hint?
> 
> The original COSMO REA6 has this SRS:
> 
> ```bash
> $ gdalinfo T_2M.2D.201501.grb
> Driver: GRIB/GRIdded Binary (.grb)
> Files: ../T_2M.2D.201501.grb
> Size is 848, 824
> Coordinate System is:
> GEOGCS["Coordinate System imported from GRIB file",
>     DATUM["unknown",
>         SPHEROID["Sphere",6367470,0]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433]]
> Origin = (-28.430500590318776,21.890500607533415)
> Pixel Size = (0.055001180637544,-0.055001215066829)
> Corner Coordinates:
> Upper Left  ( -28.4305006,  21.8905006) ( 28d25'49.80"W, 21d53'25.80"N)
> Lower Left  ( -28.4305006, -23.4305006) ( 28d25'49.80"W, 23d25'49.80"S)
> Upper Right (  18.2105006,  21.8905006) ( 18d12'37.80"E, 21d53'25.80"N)
> Lower Right (  18.2105006, -23.4305006) ( 18d12'37.80"E, 23d25'49.80"S)
> Center      (  -5.1100000,  -0.7700000) (  5d 6'36.00"W,  0d46'12.00"S)
> 
> ```
> 
> Since the original file is a bit heavy (~700Mb), I've extracted the
> first GeoTif using:
> 
> ```bash
> $ gdal_translate -b 1 -a_nodata 9999 -co "ZLEVEL=9" -co
> "COMPRESS=DEFLATE" -co "PREDICTOR=1" -co "TILED=YES" -co "TFW=YES"
> T_2M.2D.201501.grb T_2m_201501_01.tif
> ```
> 
> The geotif is available here:
> 
> https://gist.github.com/zarch/ad2274aec9746dc6267836b39f72d145
> 
> I put it into a gist to make the whole procedure reproducible.
> 
> The geotiff preserved the SRS:
> 
> ```bash
> $ gdalsrsinfo T_2m_201501_01.tif
> 
> PROJ.4 : '+proj=longlat +a=6367470 +b=6367470 +no_defs '
> 
> OGC WKT :
> GEOGCS["Coordinate System imported from GRIB file",
>     DATUM["unknown",
>         SPHEROID["Sphere",6367470,0]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433]]
> ```
> 
> I created a Location using the tif file:
> 
> ```bash
> $ grass73 -c T_2m_201501_01.tif cosmo_rea6_original
> ```
> 
> The location seems to have the same projection too:
> 
> ```bash
> $ g.proj -j
> +proj=longlat
> +no_defs
> +a=6367470
> +b=6367470
> 
> $ g.proj -w
> GEOGCS["unnamed",
>     DATUM["unknown",
>         SPHEROID["unnamed",6367470,"inf"]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433]]
> ```
> 
> Therefore I imported the raster file using r.in.gdal:
> 
> ```bash
> $ r.in.gdal input=T_2m_201501_01.tif output=T_2m_201501_01 memory=2047
> title="COSMO REA6"
> ```
> 
> Then I created a new location using the EPSG:3035
> 
> ```bash
> $ grass73 -c EPSG:3035 cosmo_rea6_3035
> ```
> 
> And I'm imported the raster map using: r.proj
> 
> ```bash
> $ g.region `r.proj location=cosmo_rea6_original mapset=PERMANENT
> input=T_2m_201501_01 -g` -p
> Input map <T_2m_201501_01 at PERMANENT> in location 
> <cosmo_rea6_original>
> :
> projection: 99 (ETRS89 / LAEA Europe)
> zone:       0
> datum:      etrs89
> ellipsoid:  grs80
> north:      869787.52797999
> south:      -4562269.74067828
> west:       -515538.58457224
> east:       5379827.69647935
> nsres:      6592.30251051
> ewres:      6952.0828786
> rows:       824
> cols:       848
> cells:      698752
> 
> $ r.proj location=cosmo_rea6_original mapset=PERMANENT
> input=T_2m_201501_01 method=nearest memory=2047
> ```
> 
> Now I downloaded the shape file with the statistical units:
> 
> ```bash
> $ wget
> http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2013_01M_SH.zip
> $ unzip NUTS_2013_01M_SH.zip
> ```
> 
> Now I created a new location using the shape file to define the location
> projection:
> 
> ```bash
> $ grass73 -c NUTS_2013_01M_SH/data/NUTS_RG_01M_2013.shp nuts_etrs89
> ```
> 
> I imported the file with:
> 
> ```bash
> $ v.in.ogr input=NUTS_2013_01M_SH/data/NUTS_RG_01M_2013.shp
> output=nusts_2013 snap=1e-05
> ```
> 
> and I reprojected the vector file into the cosmo_rea6_3035 location.
> 
> ```bash
> v.proj location=nuts_etrs89 mapset=PERMANENT input=nusts_2013
> ```
> 
> Now If I display the COSMO REA6 raster data set and the NUST3 vector
> data set I see the attached figure.
> 
> Do you know the right way to properly re-project the COSMO REA6 data set
> into a EPSG3035 location?
> 
> Best regards
> 
> Pietro
> 
> 
> _______________________________________________
> grass-user mailing list

> grass-user at .osgeo

> https://lists.osgeo.org/mailman/listinfo/grass-user
> 
> not_matching.png (86K)
> <http://osgeo-org.1560.x6.nabble.com/attachment/5341338/0/not_matching.png>

not really an answer, but the data documentation says:

"Output format: 	GRIB 1, rotated spherical"

a quick web search shows some hints to handle such data.



-----
best regards
Helmut
--
Sent from: http://osgeo-org.1560.x6.nabble.com/Grass-Users-f3884509.html