[R-sig-Geo] Problem with projection in NetCDF file

[Oscar Perpiñán Lamigueiro]

Hello,

Sorry in advance for this long email.  

I am trying to use the NetCDF files available at the
Meteogalicia-THREDDS server
(http://www.meteogalicia.es/web/modelos/threddsIndex.action). I am
having problems with the projection parameters (see figure attached for
a quick understanding of the problem).

For example, I download and open the file
wrf_arw_det_history_d01_20130214_0000.nc4 (~85Mb):

    library(raster) ## 2.1.6
    library(rgdal) ## 0.8.1

    download.file('http://mandeo.meteogalicia.es/thredds/fileServer/wrf_2d_36km/fmrc/files/20130214/wrf_arw_det_history_d01_20130214_0000.nc4',
                  destfile='~/temp/wrf_arw_det_history_d01_20130214_0000.nc4')

    ## Temperature values
    temps <- brick('~/temp/wrf_arw_det_history_d01_20130214_0000.nc4', varname='temp')
    ## This file uses the LCC projection
    projection(temps)
    ## "+proj=lcc +lon_0=-14.1000003814697 +lat_0=24.2280006408691
    ## +x_0=2182.62935 +y_0=-269.65597 +lat_1=43 +lat_2=43 +ellps=WGS84"

    ## Latitude
    lats <- raster('~/temp/wrf_arw_det_history_d01_20130214_0000.nc4',
                   varname='lat')
    ## However the latitude layer does not inherit that projection...
    projection(lats)
    ## "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
    ## Therefore, a correction is needed
    projection(lats) <- projection(temps)

    ## Longitude
    lons <- raster('~/temp/wrf_arw_det_history_d01_20130214_0000.nc4', 
                    varname='lon')
    projection(lons) <- projection(temps)

I detect the problem with the projection parameters when overlaying the
administrative boundaries:

    download.file('http://gadm.org/data/rda/ESP_adm0.RData',
              destfile='~/temp/ESP_adm0.RData')
    load('~/temp/ESP_adm0.RData') ## a new object named gadm
    projection(gadm)
    ## "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

    projLCC <- projection(temps)
    spainLCC <- spTransform(gadm, CRS(projLCC))

    library(rasterVis)
    levelplot(temps, layers=10) +
      layer(sp.polygons(spainLCC, lwd=.4))+
      contourplot(lons) + contourplot(lats)

Then, after inspecting the metadata I found that x and y are using km
units. Since this information has not been transferred to the proj4
string, I modify the projection:

    projLCC <- paste(projection(temps), '+units=km')
    projection(temps) <- projLCC
    projection(lons) <- projLCC
    projection(lats) <- projLCC
    spainLCC <- spTransform(gadm, CRS(projLCC))

However, the result is still not correct (see temps.jpg):

    levelplot(temps, layers=10) +
      layer(sp.polygons(spainLCC, lwd=.4))+
      contourplot(lons) + contourplot(lats)

Could it be that the parameters of the proj4 string are incorrect? Maybe
am I missing something?  I have tested this same file with the Godiva2
viewer and it is correctly displayed:
http://mandeo.meteogalicia.es/thredds/godiva2/godiva2.html?menu=&layer=temp&elevation=0&time=2013-02-22T14:29:45Z&scale=267.15,313.15&bbox=-49.182593,13.500592,18.788996,66.603396&server=http://mandeo.meteogalicia.es/thredds/wms/wrf_2d_36km/fmrc/files/20130214/wrf_arw_det_history_d01_20130214_0000.nc4.

Thanks in advance.

Best,

Oscar.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: temps.jpg
Type: image/jpeg
Size: 19241 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130222/cf43dbea/attachment.jpg>
-------------- next part --------------
-- 
Oscar Perpi??n Lamigueiro
Grupo de Sistemas Fotovoltaicos (IES-UPM)
Dpto. Ingenier?a El?ctrica (EUITI-UPM)
URL: http://procomun.wordpress.com
Twitter: @oscarperpinan