[R-sig-Geo] Small bug in raster reading a NetCDF?

Pascal Oettli kridox at ymail.com
Mon May 27 08:53:04 CEST 2013


Hello,

Extent and length of latitudes are correct.

In NCEP/NCAR reanalysis, coordinates are the centers of grid cells. 
Thus, for North Pole, the grid extents (+---+), centers (o) and Earth 
"limits" (____) are:


91.25 +-------+-------+-------+-------+
       |               |               | 1.25
       |               |               |
__90__+_______o_______|_______o_______|______
       |               |               |
       |               |               | 1.25
88.75 +-------+-------+-------+-------+
       |               |               | 1.25
       |               |               |
  87.5 +       o       +       o       +
       |               |               |
       |               |               | 1.25
86.25 +-------+-------+-------+-------+
     -1.25     0      1.25    2.5     3.75


180 = 73 x 2.5 - (2 x 1.25)
                       ^
             |91.25-90| + |-91.25+90|

Regards,
Pascal


On 05/26/2013 07:13 AM, Jonathan Greenberg wrote:
> Curiouser and curiouser -- I feel like I'm missing some obvious point here
> -- the NCEP grid can't fit into its specified extent, pixel size, and
> number of pixels -- the y-axis seems to have one pixel too many (e.g. it
> has an nrow of 73, but 73 x 2.5 degrees =  182.5 degrees, which can't be
> true (the specified pixel size from NOAA,
> http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.surface.html).
>   Does anyone have any idea whether this is an "extra row", or if the pixel
> sizes are not set properly, or some other point I'm failing to understand?
>
> --j
>
>
> On Sat, May 25, 2013 at 4:52 PM, Greenberg, Jonathan <jgrn at illinois.edu>wrote:
>
>> Robert et al:
>>
>> I think there may be a small bug in raster reading an NCEP NetCDF file.  I
>> think the NetCDF file is supplying the upper left coordinates of each cell,
>> but raster is assuming these coordinates are the center of the cell:
>>
>> # Given:
>> require("ncdf")
>> require("raster")
>> # Sample global netCDF from NCEP Reanalysis, about a 30mb download:
>> download.file("
>> ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface/air.sig995.1948.nc
>> ","air.sig995.1948.nc",mode="wb")
>> # Create a brick:
>> myncdf <- brick("air.sig995.1948.nc")
>> myncdf
>> # Notice the extent has impossible coordinates:
>> class       : RasterBrick
>> dimensions  : 73, 144, 10512, 1464  (nrow, ncol, ncell, nlayers)
>> resolution  : 2.5, 2.5  (x, y)
>> extent      : -1.25, 358.75, -91.25, 91.25  (xmin, xmax, ymin, ymax)
>> coord. ref. : +proj=longlat +datum=WGS84
>> data source : /Users/jgrn307/air.sig995.1948.nc
>> names       : X1948.01.03.00.00.00, X1948.01.03.06.00.00,
>> X1948.01.03.12.00.00, X1948.01.03.18.00.00, X1948.01.04.00.00.00,
>> X1948.01.04.06.00.00, X1948.01.04.12.00.00, X1948.01.04.18.00.00,
>> X1948.01.05.00.00.00, X1948.01.05.06.00.00, X1948.01.05.12.00.00,
>> X1948.01.05.18.00.00, X1948.01.06.00.00.00, X1948.01.06.06.00.00,
>> X1948.01.06.12.00.00, ...
>> Date/time   : 1948-01-03 00:00:00, 1949-01-02 18:00:00 (min, max)
>> varname     : air
>>
>> # Opening the NetCDF directly:
>> myncdf_con <- open.ncdf("air.sig995.1948.nc")
>> get.var.ncdf(myncdf_con,"lon")
>>
>> # Output:
>>    [1]   0.0   2.5   5.0   7.5  10.0  12.5  15.0  17.5  20.0  22.5  25.0
>>   27.5  30.0  32.5  35.0  37.5  40.0  42.5  45.0  47.5  50.0  52.5
>>   [23]  55.0  57.5  60.0  62.5  65.0  67.5  70.0  72.5  75.0  77.5  80.0
>>   82.5  85.0  87.5  90.0  92.5  95.0  97.5 100.0 102.5 105.0 107.5
>>   [45] 110.0 112.5 115.0 117.5 120.0 122.5 125.0 127.5 130.0 132.5 135.0
>> 137.5 140.0 142.5 145.0 147.5 150.0 152.5 155.0 157.5 160.0 162.5
>>   [67] 165.0 167.5 170.0 172.5 175.0 177.5 180.0 182.5 185.0 187.5 190.0
>> 192.5 195.0 197.5 200.0 202.5 205.0 207.5 210.0 212.5 215.0 217.5
>>   [89] 220.0 222.5 225.0 227.5 230.0 232.5 235.0 237.5 240.0 242.5 245.0
>> 247.5 250.0 252.5 255.0 257.5 260.0 262.5 265.0 267.5 270.0 272.5
>> [111] 275.0 277.5 280.0 282.5 285.0 287.5 290.0 292.5 295.0 297.5 300.0
>> 302.5 305.0 307.5 310.0 312.5 315.0 317.5 320.0 322.5 325.0 327.5
>> [133] 330.0 332.5 335.0 337.5 340.0 342.5 345.0 347.5 350.0 352.5 355.0
>> 357.5
>>
>> get.var.ncdf(myncdf_con,"lat")
>> # Output:
>>   [1]  90.0  87.5  85.0  82.5  80.0  77.5  75.0  72.5  70.0  67.5  65.0
>>   62.5  60.0  57.5  55.0  52.5  50.0  47.5  45.0  42.5  40.0  37.5
>> [23]  35.0  32.5  30.0  27.5  25.0  22.5  20.0  17.5  15.0  12.5  10.0
>> 7.5   5.0   2.5   0.0  -2.5  -5.0  -7.5 -10.0 -12.5 -15.0 -17.5
>> [45] -20.0 -22.5 -25.0 -27.5 -30.0 -32.5 -35.0 -37.5 -40.0 -42.5 -45.0
>> -47.5 -50.0 -52.5 -55.0 -57.5 -60.0 -62.5 -65.0 -67.5 -70.0 -72.5
>> [67] -75.0 -77.5 -80.0 -82.5 -85.0 -87.5 -90.0
>>
>> # Don't forget to close the connection:
>> close.ncdf(myncdf_con)
>>
>> Thoughts?
>>
>> --j
>>
>> --
>> Jonathan A. Greenberg, PhD
>> Assistant Professor
>> Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
>> Department of Geography and Geographic Information Science
>> University of Illinois at Urbana-Champaign
>> 607 South Mathews Avenue, MC 150
>> Urbana, IL 61801
>> Phone: 217-300-1924
>> http://www.geog.illinois.edu/~jgrn/
>> AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
>>
>>          [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>



More information about the R-sig-Geo mailing list