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

[Pietro Zambelli]

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not_matching.png
Type: image/png
Size: 64842 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20171106/c3fc0d9a/attachment-0001.png>