I am working with WRF climate data (4km resolution for North America, available here: https://rda.ucar.edu/datasets/ds612.0/index.html#!description). The data format is netcdf. I would like to read the data in as a raster and project the data into lat/lon coordinates using R. I would like to do this so I can extract spatial subsets of the data and use raster algebra. However, I am having trouble projecting the data.
There has been a lot written online about WRF projections and how tricky they are to work with (see this and this).
To start with I've been using the file of constants (lat, lon, terrain height, soil type, etc). It is available here from the original source: https://rda.ucar.edu/datasets/ds612.0/index.html#cgi-bin/datasets/getWebList?dsnum=612.0&action=customize&disp=&gindex=1, it is called RALconus_4km_wrf_constants.nc. I've made a copy of the file available here. This file contains a field 'HGT' (terrain height) with dimensions [west_east,south_north,Time] as well as fields 'XLAT' and 'XLONG' with dimensions [west_east,south_north]. The latitude and longitude values are unique for each grid cell. From print(nc):
Time Size:1 *** is unlimited ***
units: hours since 1901-01-01 00:00:00
calendar: standard
long_name: Time
description: Time
south_north Size:1015
west_east Size:1359
I know from the WRF documentation (here), that the projection is lambert conformal conic, and the netcdf file provides a few additional details:
110 global attributes:
TITLE: OUTPUT FROM WRF V3.4.1 MODEL
START_DATE: 2001-02-28_00:00:00
SIMULATION_START_DATE: 2000-10-01_00:00:00
WEST-EAST_GRID_DIMENSION: 1360
SOUTH-NORTH_GRID_DIMENSION: 1016
BOTTOM-TOP_GRID_DIMENSION: 51
DX: 4000
DY: 4000
GRIDTYPE: C
DIFF_OPT: 1
KM_OPT: 4
DAMP_OPT: 3
DAMPCOEF: 0.200000002980232
KHDIF: 0
KVDIF: 0
MP_PHYSICS: 28
RA_LW_PHYSICS: 4
RA_SW_PHYSICS: 4
SF_SFCLAY_PHYSICS: 11
SF_SURFACE_PHYSICS: 4
BL_PBL_PHYSICS: 1
CU_PHYSICS: 0
SURFACE_INPUT_SOURCE: 1
SST_UPDATE: 1
GRID_FDDA: 2
GFDDA_INTERVAL_M: 360
GFDDA_END_H: 999999
GRID_SFDDA: 0
SGFDDA_INTERVAL_M: 0
SGFDDA_END_H: 0
HYPSOMETRIC_OPT: 2
SF_URBAN_PHYSICS: 0
SHCU_PHYSICS: 0
MFSHCONV: 0
FEEDBACK: 0
SMOOTH_OPTION: 2
SWRAD_SCAT: 1
W_DAMPING: 1
MOIST_ADV_OPT: 1
SCALAR_ADV_OPT: 1
TKE_ADV_OPT: 1
DIFF_6TH_OPT: 0
DIFF_6TH_FACTOR: 0.119999997317791
FGDT: 0
GUV: 4.99999987368938e-05
GT: 4.99999987368938e-05
GPH: 4.99999987368938e-05
IF_RAMPING: 1
DTRAMP_MIN: 60
OBS_NUDGE_OPT: 0
BUCKET_MM: 100
BUCKET_J: 1e+09
PREC_ACC_DT: 60
OMLCALL: 0
ISFTCFLX: 0
ISHALLOW: 0
OPT_SFC: 1
DVEG: 4
OPT_CRS: 1
OPT_BTR: 2
OPT_RUN: 1
OPT_FRZ: 1
OPT_INF: 1
OPT_RAD: 3
OPT_ALB: 2
OPT_SNF: 4
OPT_TBOT: 1
OPT_STC: 1
DFI_OPT: 0
WEST-EAST_PATCH_START_UNSTAG: 1
WEST-EAST_PATCH_END_UNSTAG: 1359
WEST-EAST_PATCH_START_STAG: 1
WEST-EAST_PATCH_END_STAG: 1360
SOUTH-NORTH_PATCH_START_UNSTAG: 1
SOUTH-NORTH_PATCH_END_UNSTAG: 1015
SOUTH-NORTH_PATCH_START_STAG: 1
SOUTH-NORTH_PATCH_END_STAG: 1016
BOTTOM-TOP_PATCH_START_UNSTAG: 1
BOTTOM-TOP_PATCH_END_UNSTAG: 50
BOTTOM-TOP_PATCH_START_STAG: 1
BOTTOM-TOP_PATCH_END_STAG: 51
GRID_ID: 1
PARENT_ID: 0
I_PARENT_START: 1
J_PARENT_START: 1
PARENT_GRID_RATIO: 1
DT: 15
CEN_LAT: 39.7000122070312
CEN_LON: -98
TRUELAT1: 28
TRUELAT2: 50
MOAD_CEN_LAT: 39.7000122070312
STAND_LON: -98
POLE_LAT: 90
POLE_LON: 0
GMT: 0
JULYR: 2001
JULDAY: 59
MAP_PROJ: 1
MMINLU: MODIFIED_IGBP_MODIS_NOAH
NUM_LAND_CAT: 20
ISWATER: 17
ISLAKE: -1
ISICE: 15
ISURBAN: 13
ISOILWATER: 14
Conventions: CF-1.6
history: Tue Apr 4 15:29:54 2017: ncks -A mylatlon4.nc myclean.nc
history_of_appended_files: Tue Apr 4 15:29:54 2017: Appended file mylatlon4.nc had no "history" attribute
NCO: "4.5.5"
I am starting by trying to project the HGT field. The raw data from the HGT field from the netcdf file looks like this.
In my attempts to properly project the raster, I have looked at this and other posts, but I have not been successful yet. I have tried the following code and variants of it:
wrfdir <- '/Volumes/Web/abby/DATA/WRF/'
metafn <- paste0(wrfdir, 'RALconus4km_wrf_constants.nc')
rr <- raster(metafn,varname='HGT')
projection(rr) <- '+proj=lcc +lat_1=28.0 +lat_2=50.0 +lat_0=39.7000122070312 +lon_0=-98.0 +ellps=sphere +a=6370000 +b=6370000 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
# then project to lat lon
rrLL = projectRaster(rr, crs="+init=epsg:4326")
plot(rrLL[[1]])
This results in the following map, in which the lats and lons are obviously incorrect: first projection attempt
I have also tried the following, based on this post:
crs(rr) <- '+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs'
rrLL <- projectRaster((rr),crs='+proj=lcc +lat_1=28.0 +lat_2=50.0 +lat_0=39.7000122070312 +lon_0=-98.0 +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs')
Which results in the following error message:
Error in if (nr != x@nrows | nc != x@ncols) { : missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
634 projected point(s) not finite
2: In rgdal::rawTransform(projection(obj), crs, nrow(xy), xy[, 1], :
4 projected point(s) not finite
Any thoughts on how to do this would be much appreciated!
