[R-sig-Geo] problems drawing a world map on a levelplot

Tom Roche Tom_Roche at pobox.com
Wed Dec 19 19:52:38 CET 2012


summary: How to best/easiest range-shift lon-lat data?

details:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017011.html
> I happened to write [the following code, omitted] yesterday, which
> switches the parts of wrld_simpl from < 0 longitude to the east
> (Pacific view).

<self-administered dopeslap/> Of course! the problem is my input!
see starred line:

R session
> > global.proj <- CRS('+proj=longlat +ellps=WGS84')
> > global.raster <- brick(in.fp, varname=data.var.name, crs=global.proj)
> > global.raster
> class       : RasterBrick 
> dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
> resolution  : 2.5, 1.894737  (x, y)
* extent      : -1.25, 358.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

which is wrong:

- latitudes outside the range [-90, +90]?

- (as you note) 0-based longitudes are likely troublesome (I'm pretty
  sure the other data with which this must integrate is [-180, +180])

- ... and these longitudes aren't even 0-based!

I didn't even check :-( since the data's source is much more eminent
than a mere student like myself. But now I note

$ ncdump -h ./2008N2O_restart.nc | head
> netcdf \2008N2O_restart {
> dimensions:
>         lev = 56 ;
>         lat = 96 ;
>         lon = 144 ;
>         ilev = 57 ;
> variables:
>         double N2O(lev, lat, lon) ;
>                 N2O:long_name = "N2O" ;
>                 N2O:units = "ppmV" ;

$ ncdump -v lat ./2008N2O_restart.nc | tail -n 28
>  lat = -90, -88.1052631578947, -86.2105263157895, -84.3157894736842, 
>     -82.4210526315789, -80.5263157894737, -78.6315789473684, 
>     -76.7368421052632, -74.8421052631579, -72.9473684210526, 
>     -71.0526315789474, -69.1578947368421, -67.2631578947368, 
>     -65.3684210526316, -63.4736842105263, -61.5789473684211, 
>     -59.6842105263158, -57.7894736842105, -55.8947368421053, -54, 
>     -52.1052631578947, -50.2105263157895, -48.3157894736842, 
>     -46.4210526315789, -44.5263157894737, -42.6315789473684, 
>     -40.7368421052632, -38.8421052631579, -36.9473684210526, 
>     -35.0526315789474, -33.1578947368421, -31.2631578947368, 
>     -29.3684210526316, -27.4736842105263, -25.5789473684211, 
>     -23.6842105263158, -21.7894736842105, -19.8947368421053, -18, 
>     -16.1052631578947, -14.2105263157895, -12.3157894736842, 
>     -10.4210526315789, -8.52631578947369, -6.63157894736842, 
>     -4.73684210526316, -2.84210526315789, -0.947368421052634, 
>     0.947368421052634, 2.84210526315788, 4.73684210526315,
>     6.63157894736842, 8.52631578947369, 10.4210526315789,
>     12.3157894736842, 14.2105263157895, 16.1052631578947, 18,
>     19.8947368421053, 21.7894736842105, 23.6842105263158,
>     25.5789473684211, 27.4736842105263, 29.3684210526316,
>     31.2631578947368, 33.1578947368421, 35.0526315789474,
>     36.9473684210526, 38.8421052631579, 40.7368421052632,
>     42.6315789473684, 44.5263157894737, 46.4210526315789,
>     48.3157894736842, 50.2105263157895, 52.1052631578947, 54,
>     55.8947368421053, 57.7894736842105, 59.6842105263158,
>     61.578947368421, 63.4736842105263, 65.3684210526316,
>     67.2631578947368, 69.1578947368421, 71.0526315789474,
>     72.9473684210526, 74.8421052631579, 76.7368421052632,
>     78.6315789473684, 80.5263157894737, 82.4210526315789,
>     84.3157894736842, 86.2105263157895, 88.1052631578947, 90 ; }

I'm not sure why raster is reporting y range=[-90.94737, 90.94737],
except that 2 * 0.94737 = 1.89474 == the latitude resolution.

$ ncdump -v lon ./2008N2O_restart.nc | tail -n 15
>  lon = 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30,
>     32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60, 62.5,
>     65, 67.5, 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95,
>     97.5, 100, 102.5, 105, 107.5, 110, 112.5, 115, 117.5, 120, 122.5,
>     125, 127.5, 130, 132.5, 135, 137.5, 140, 142.5, 145, 147.5, 150,
>     152.5, 155, 157.5, 160, 162.5, 165, 167.5, 170, 172.5, 175, 177.5,
>     180, 182.5, 185, 187.5, 190, 192.5, 195, 197.5, 200, 202.5, 205,
>     207.5, 210, 212.5, 215, 217.5, 220, 222.5, 225, 227.5, 230, 232.5,
>     235, 237.5, 240, 242.5, 245, 247.5, 250, 252.5, 255, 257.5, 260,
>     262.5, 265, 267.5, 270, 272.5, 275, 277.5, 280, 282.5, 285, 287.5,
>     290, 292.5, 295, 297.5, 300, 302.5, 305, 307.5, 310, 312.5, 315,
>     317.5, 320, 322.5, 325, 327.5, 330, 332.5, 335, 337.5, 340, 342.5,
>     345, 347.5, 350, 352.5, 355, 357.5 ; }

Again, it seems raster is reporting range(x)=[-1.25, 358.75] only
because resolution(x)=2.5°. I'll defer to its superior wisdom.

Given my need to integrate this data with other more conventional,
longitudes=[-180, +180] datasets, it seems better to fix this data than
to fix wrld_simpl. So I'm wondering, how best to "shift" my data from
longitudes=[0, 360] to longitudes=[-180, +180]?

TIA, Tom Roche <Tom_Roche at pobox.com>



More information about the R-sig-Geo mailing list