3

I have a raster, obtained from a netcdf which is in (Lambert Conic Conformal projection):

    library(meteoForecast)
    wrf_temporary <- getRaster("temp", day = Sys.Date(), frames = 'complete', resolution = 36, service = "meteogalicia")
    wrf_temporary

extent      : -18, 4230, -18, 3726  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=43 +lat_2=43 +lat_0=34.82300186157227 +lon_0=-14.10000038146973 +x_0=536402.34 +y_0=-18558.61 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs

Now I want to transform that wrf_temporary raster to "+proj=longlat +datum=WGS84" (lat long degree). What to do? I want something like this:

    mfExtent('meteogalicia', resolution = 36)

class       : Extent 
xmin        : -49.18259 
xmax        : 18.789 
ymin        : 24.03791 
ymax        : 56.06608

Already tried a lot of options but none of them give the right results...

Ricardo Faria
  • 111
  • 9

1 Answers1

4

This should work:

> ll = projectRaster(wrf_temporary,crs="+init=epsg:4326")
> plot(ll[[1]])

and produces this:

enter image description here

which looks right at first glance but the Greenwich meridian (longitude=0) doesn't go through London.

The problem doesn't occur for resolution set to the other values, 4 or 12. When called with resolution=36, you get a raster with a larger extent than the others, but with an identical projection string. This can't be right. For example, the (0,0) coordinate on the following plots should be the same point on the earth because they claim to be the same projection, but they arent.

enter image description here

I think the resolution=12 ones are correct. Here's that raster transformed to lat-long with an EU vector coverage overlaid:

enter image description here

which lines up perfectly.

So I would say its a bug - either getRaster is guessing the projection and getting it wrong, or not applying a change in projection after cropping, or whatever service it uses is lying about the projection.

getRaster is guessing. The projection info is in the downloaded NetCDF file. In another format. For the resolution=36 file the projection info section is this:

int Lambert_Conformal ;
    Lambert_Conformal:standard_parallel = 43., 43. ;
    Lambert_Conformal:latitude_of_projection_origin = 24.2280006408691 ;
    Lambert_Conformal:false_easting = 2182.62935 ;
    Lambert_Conformal:false_northing = -269.65597 ;
    Lambert_Conformal:grid_mapping_name = "lambert_conformal_conic" ;
    Lambert_Conformal:longitude_of_central_meridian = -14.1000003814697 ;
    Lambert_Conformal:_CoordinateTransformType = "Projection" ;
    Lambert_Conformal:_CoordinateAxisTypes = "GeoX GeoY" ;

which is a bit different to the resolution=12 projection used by the R package. So make a compliant one:

p = "+proj=lcc +lat_1=43 +lat_2=43 +lat_0=24.2280006408691 +lon_0=-14.1000003814697 +x_0=2182629.35 +y_0=-269655.97 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs"
 wrf <- getRaster('temp', day = testDay, resolution=36)
 projection(wrf) = p

Then my test with EU overlay....

enter image description here

Note this time I've reprojected the EU. Doing projectRaster to latlong works too:

> wrfLL = projectRaster(wrf, crs="+init=epsg:4326")
> plot(wrfLL[[1]])
> abline(v=0)

enter image description here

with the Greenwich meridian in its rightful place.

Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • Thanks a lot for your comment, I already tried that way, but there is a problem, that projection is not correct it is not representing lat long degree, to give you one idea of how bad it is, try to run the following code: `image(ll, layers = 1)` `plot(getMap(resolution = "high"), add = T)` Thanks, Ricardo Faria – Ricardo Faria Apr 26 '16 at 15:47
  • Ah ha. Greenwich meridian (x=0) doesn't go through London... Should have noticed... Hmmm. – Spacedman Apr 26 '16 at 15:49
  • The problem is that I need domain 03 which is the resolution=36, I need the data on Madeira and Canary islands. The resolution=12 don't cover the area I need. – Ricardo Faria Apr 26 '16 at 16:32
  • Have done some detective work for you! I see your bug report now too! – Spacedman Apr 26 '16 at 16:44
  • My fault. I was using the same projection string for the three resolutions. It is [fixed now](https://github.com/oscarperpinan/meteoForecast/commit/69f9439723fca8759ff2773fdb1a709fbbc5f834). Thanks a lot for your detailed answer. – Oscar Perpiñán Apr 26 '16 at 17:53
  • I want to thank you both for helping me with this situation. :D – Ricardo Faria Apr 26 '16 at 19:06