title: Trouble projecting WRF data in R
body: <p>I am working with WRF climate data (4km resolution for North America, available here: <a href="https://rda.ucar.edu/datasets/ds612.0/index.html#!description" rel="nofollow noreferrer">https://rda.ucar.edu/datasets/ds612.0/index.html#!description</a>).  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. </p>
  
  <p>There has been a lot written online about WRF projections and how tricky they are to work with (see <a href="https://fabienmaussion.info/2018/01/06/wrf-projection/" rel="nofollow noreferrer">this</a> and <a href="https://journals.ametsoc.org/doi/full/10.1175/MWR-D-12-00351.1" rel="nofollow noreferrer">this</a>).  </p>
  
  <p>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: <a href="https://rda.ucar.edu/datasets/ds612.0/index.html#cgi-bin/datasets/getWebList?dsnum=612.0&amp;action=customize&amp;disp=&amp;gindex=1" rel="nofollow noreferrer">https://rda.ucar.edu/datasets/ds612.0/index.html#cgi-bin/datasets/getWebList?dsnum=612.0&amp;action=customize&amp;disp=&amp;gindex=1</a>, it is called RALconus_4km_wrf_constants.nc.   I've made a copy of the file available <a href="https://drive.google.com/file/d/1ipHVP9_uF3Y-9IfIq52OdsDgJroPgLW9/view?usp=sharing" rel="nofollow noreferrer">here</a>. 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):</p>
  
  <pre><code>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
  </code></pre>
  
  <p>I know from the WRF documentation (<a href="http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm" rel="nofollow noreferrer">here</a>), that the projection is lambert conformal conic, and the netcdf file provides a few additional details:</p>
  
  <pre><code>    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"
  </code></pre>
  
  <p>I am starting by trying to project the HGT field. The raw data from the HGT field from the netcdf file looks like <a href="https://i.stack.imgur.com/FihpU.jpg" rel="nofollow noreferrer">this</a>.</p>
  
  <p><a href="https://i.stack.imgur.com/FihpU.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/FihpU.jpg" alt="enter image description here"></a></p>
  
  <p>In my attempts to properly project the raster, I have looked at <a href="https://stackoverflow.com/questions/36868506/how-to-change-a-lambert-conic-conformal-raster-projection-to-latlon-degree-r/36869273#36869273">this</a> and other posts, but I have not been successful yet.  I have tried the following code and variants of it:</p>
  
  <pre><code>  wrfdir &lt;- '/Volumes/Web/abby/DATA/WRF/'
    metafn &lt;- paste0(wrfdir, 'RALconus4km_wrf_constants.nc')
    rr &lt;- raster(metafn,varname='HGT')
    projection(rr) &lt;- '+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]])
  </code></pre>
  
  <p>This results in the following map, in which the lats and lons are obviously incorrect:
  <a href="https://i.stack.imgur.com/YjQAs.jpg" rel="nofollow noreferrer">first projection attempt</a></p>
  
  <p><a href="https://i.stack.imgur.com/YjQAs.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YjQAs.jpg" alt="enter image description here"></a></p>
  
  <p>I have also tried the following, based on this <a href="http://r-sig-geo.2731867.n2.nabble.com/adapting-spatial-points-and-wrld-smpl-to-a-reference-system-implicit-in-a-nc-file-tt7589570.html" rel="nofollow noreferrer">post</a>:</p>
  
  <pre><code> crs(rr) &lt;- '+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs' 
   rrLL &lt;- 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')
  </code></pre>
  
  <p>Which results in the following error message:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>Any thoughts on how to do this would be much appreciated!</p>