{"action":"opened","issue":{"url":"https://api.github.com/repos/ran-ran/MyNotes/issues/5","repository_url":"https://api.github.com/repos/ran-ran/MyNotes","labels_url":"https://api.github.com/repos/ran-ran/MyNotes/issues/5/labels{/name}","comments_url":"https://api.github.com/repos/ran-ran/MyNotes/issues/5/comments","events_url":"https://api.github.com/repos/ran-ran/MyNotes/issues/5/events","html_url":"https://github.com/ran-ran/MyNotes/issues/5","id":276932575,"number":5,"title":"读取grib2文件","user":{"login":"ran-ran","id":13953953,"avatar_url":"https://avatars0.githubusercontent.com/u/13953953?v=4","gravatar_id":"","url":"https://api.github.com/users/ran-ran","html_url":"https://github.com/ran-ran","followers_url":"https://api.github.com/users/ran-ran/followers","following_url":"https://api.github.com/users/ran-ran/following{/other_user}","gists_url":"https://api.github.com/users/ran-ran/gists{/gist_id}","starred_url":"https://api.github.com/users/ran-ran/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/ran-ran/subscriptions","organizations_url":"https://api.github.com/users/ran-ran/orgs","repos_url":"https://api.github.com/users/ran-ran/repos","events_url":"https://api.github.com/users/ran-ran/events{/privacy}","received_events_url":"https://api.github.com/users/ran-ran/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":0,"created_at":"2017-11-27T08:14:25Z","updated_at":"2017-11-27T08:14:25Z","closed_at":null,"author_association":"OWNER","body":"这个笔记整理了如何读取grib2格式文件。\r\n下面代码中的cmd字符串也可以在命令行工具中直接输入，效果相同。\r\n\r\n首先需要下载wgrib2命令行程序，这里是[链接](ftp://ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/Windows_XP)。将其中的wgrib2.exe和其他所有dll文件下载到某个本地路径下，并且将这个路径添加到系统环境变量中（可能需要重启生效）。\r\n\r\n完成上面步骤后，可以在windows打开cmd，输入wgrib2，看看能否运行。\r\n\r\n下面的例子中构造的cmd命令都是在R中用system函数执行的，也可以在windows的cmd中直接输入。\r\n也可以使用windows命令行的升级版powershell。\r\n\r\n下面的例子可以参考一下，特别是输入到csv文件。\r\n\r\n\r\n#显示文件清单(inventory)\r\n\r\n```{r,echo=TRUE}\r\n# 这个fname是我放文件的路径\r\nfname <- 'F:/tfiles/grib2_test/z500.ensm.jan.cfsv2.data.grb2'\r\ncmd <- paste('wgrib2 ', fname)\r\nprint(cmd)\r\ninv <- system(cmd, intern = TRUE)\r\nprint(inv[1:10])\r\n```\r\n\r\n例如上面清单中最后一行：\r\n\r\n_280:13184074:d=2009010618:HGT:500 mb:9-10 month ave fcst:_\r\n\r\n它的格式含义如下：\r\n\r\n_(msg|msg.submsg):(byte location):d=(YYYYMMDDHH):(variable):(level):(ftime):(optional attr)_\r\n\r\n+    msg = message/record number starting from 1\r\n\r\n+    submsg = submessage number starting from 1 (optional)\r\n\r\n+    byte location: byte location of the start of the message starting from 0\r\n\r\n+    YYYYMMDDHH = reference time of the grib message, YYYY = year, MM = month (01..12), DD = day, HH = hour\r\n\r\n+    level = level/layer attributes, ex. \"500 mb\", \"ground\", \"2 m above ground\"\r\n\r\n+    ftime = time attributes, ex. \"12 hour fcst\", \"anl\", \"0-6 hr ave fcst\"\r\n\r\n+    optional attr: optional attributes that are needed to define the field such as ensemble memeber, chemical type (for concentration), probability range\r\n\r\n上面的例子中没有最后的optional attr。\r\n\r\n\r\n#使用正则表达式提取一部分数据\r\n\r\n```{r,echo=TRUE,message=TRUE}\r\ncmd <- paste('wgrib2', fname, '-match \":(9-10 month ave fcst:)\"')\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\n```\r\n\r\n\r\n#使用for提取指定的数据\r\n```{r,echo=TRUE,message=TRUE}\r\ncmd <- paste('wgrib2', fname, '-for 1:3')\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\n\r\ncmd <- paste('wgrib2', fname, '-for 1:90:9')\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\n```\r\n\r\n\r\n#打印某个格点的值\r\n```{r,echo=TRUE}\r\ncmd <- paste('wgrib2', fname, '-ij 1 1 -match \":(9-10 month ave fcst:)\"')\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\n```\r\n\r\n\r\n#输出为csv文件\r\n```{r,echo=TRUE,message=TRUE,fig.height=8,fig.width=10}\r\nfout <- 'F:/tfiles/grib2_test/test.csv'\r\ncmd <- paste('wgrib2', fname, '-match \":(6-7 month ave fcst:)\"', '-match \":(d=2009010618)\"', '-csv', fout)\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\ndata <- read.csv(fout, header = F)\r\nhead(data)\r\nlibrary(lattice)\r\ncontourplot(V7~V5+V6, data = data, region = T, cuts = 20)\r\n```\r\n\r\n\r\n#输出为netcdf文件\r\n```{r,echo=TRUE,message=TRUE}\r\nfout <- 'F:/tfiles/grib2_test/test.nc'\r\ncmd <- paste('wgrib2', fname, '-match \":(6-7 month ave fcst:)\"', '-netcdf', fout)\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\nlibrary(ncdf4)\r\nnc <- nc_open(fout)\r\n```\r\n\r\n\r\n# 输出小范围的grib2文件\r\n可以用下面的方法输出一个子区域文件：\r\n\r\n_-small_grib LonW:LonE LatS:LatN file_name_\r\n\r\n```{r,echo=TRUE,fig.width=10,fig.height=8}\r\n# 首先提取区域girb2文件\r\ngrib.easia <- 'F:/tfiles/grib2_test/eaisa.grb2'\r\ncmd <- paste('wgrib2', fname, '-match \":(6-7 month ave fcst:)\" -small_grib 60:140 0:60', grib.easia)\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\n\r\n# 显示新grib2文件的经纬度范围\r\ncmd <- paste('wgrib2', grib.easia, '-domain')\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\n\r\n# 输出这个文件为csv文件\r\ncsv.easia <- 'F:/tfiles/grib2_test/easia.csv'\r\ncmd <- paste('wgrib2', grib.easia, '-csv', csv.easia)\r\ncat(cmd)\r\nsystem(cmd, intern = TRUE)\r\ndata <- read.csv(csv.easia, header = F)\r\nhead(data)\r\nlibrary(lattice)\r\nwindows()\r\ncontourplot(V7~V5+V6|V1, data = data, subset = V1 %in% c('1998-01-06 18:00:00', '1983-01-06 18:00:00', '1991-01-06 18:00:00', '1985-01-06 18:00:00', '2008-01-06 18:00:00', '2006-01-06 18:00:00'), region = T, cuts = 14)\r\n```\r\n\r\n\r\n\r\n"}}
{"actor":{"display_login":"ran-ran"}}