An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130525/add5866f/attachment.pl>
Small bug in raster reading a NetCDF?
3 messages · Jonathan Greenberg, Ping Yang, Pascal Oettli
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130525/d2b84d44/attachment.pl>
1 day later
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