Skip to content

Set up raster data in WRF Mercator projection

4 messages · David Wang, Dominik Schneider

Original
David Wang
# 
Hello,


I'd like to load a WRF (Weather Research and Forecasting) output file in NetCDF as rasters and set up correct projection and coordinates. Using land mask as an example,


landmask <- raster("geo_em_d01.nc", varname = "LANDMASK")


The projection of this particular WRF simulation is Mercator. After extensive googling, I came up with this projection string:


projection(landmask) <- CRS("+proj=merc +lat_ts=30 +lon_0=-67 +a=6370000.0 +b=6370000.0 +units=km +no_defs")


That is based on the grid information from ncdump -h geo_em_d01.nc:


                :GRIDTYPE = "C" ;
                :DX = 27000.f ;
                :DY = 27000.f ;
                :DYN_OPT = 2 ;
                :CEN_LAT = 32.f ;
                :CEN_LON = -67.f ;
                :TRUELAT1 = 30.f ;
                :TRUELAT2 = 45.f ;
                :MOAD_CEN_LAT = 32.f ;
                :STAND_LON = 290.f ;
                :POLE_LAT = 90.f ;
                :POLE_LON = 0.f ;
                :corner_lats = 6.033165f, 52.29568f, 52.29568f, 6.033165f, 6.033165f, 52.29568f, 52.29568f, 6.033165f, 5.893715f, 52.38135f, 52.38135f, 5.893715f, 5.893715f, 52.38135f, 52.38135f, 5.893715f ;
                :corner_lons = -129.8151f, -129.8151f, -4.184856f, -4.184856f, -129.9554f, -129.9554f, -4.044647f, -4.044647f, -129.8151f, -129.8151f, -4.184856f, -4.184856f, -129.9554f, -129.9554f, -4.044647f, -4.044647f ;
                :MAP_PROJ = 3 ;

However, I'm not entirely sure this projection is correct. Furthermore, I have no clue how to set extent(landmask) in km. Can anyone familiar with WRF give me a pointer or two?


Thanks,

David
Dominik Schneider
# 
did you see this thread? might help:
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
Projection_String = ('PROJCS["Sphere_Mercator",'
                             'GEOGCS["GCS_Sphere",'
                             'DATUM["D_Sphere",'
                             'SPHEROID["Sphere",6370000.0,0.0]],'
                             'PRIMEM["Greenwich",0.0],'
                             'UNIT["Degree",0.0174532925199433]],'
                             'PROJECTION["Mercator"],'
                             'PARAMETER["False_Easting",0.0],'
                             'PARAMETER["False_Northing",0.0],'
                             'PARAMETER["Central_Meridian",' +
str(central_meridian) + '],'
                             'PARAMETER["Standard_Parallel_1",' +
str(standard_parallel_1) + '],'
                             'UNIT["Meter",1.0]]')


where central_meridian = STAND_LON and standard_parallel_1=TRUELAT1

so to start I think you want (changes in *bold*):
"+proj=merc +lat_ts=30 +lon_0=*290* +a=6370000.0 +b=6370000.0 +units=*m*
 +no_defs")

if you post the data, someone might be able to help create a proper raster
for you.
On Wed, Aug 3, 2016 at 9:46 AM, David Wang <dw2116 at outlook.com> wrote:

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      David Wang
    

    

      
      #
    

  


  

      
Hello Dominik,


Thanks for your comment. Yes, I have read the thread that you pointed to prior to posting the question. It's Lambert conformal conic projection though. I'm not sure the way to compute the domain extent is the same.


I have uploaded a copy of the grid file at https://www.dropbox.com/s/ffrxqkxr4vai8hf/geo_em.d01.nc?dl=0 if anyone happens to be able to take a look.

Thanks,
David

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Dominik Schneider
    

    

      
      #
    

  


  

      
Same process, just different projection.

with this folder structure:
project_dir
- data
- scripts

the following script in 'scripts', nc file in 'data'


library(ncdf4)
library(raster)

# open the nc file for reading
ncid=nc_open('../data/geo_em.d01.nc')

#check which grid the variable of interest is on
ncatt_get(ncid,'LANDMASK')#stagger=M

# domain corner grids. see
http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm
cornerlat=ncatt_get(ncid,0)$corner_lats
cornerlong=ncatt_get(ncid,0)$corner_lons

# projected grid resolution
dx=ncatt_get(ncid,0)$DX
dy=ncatt_get(ncid,0)$DY

# the coordinates are on a sphere in longlat
proj1=CRS('+proj=longlat +towgs84=0,0,0 +a=6370000 +b=6370000 +units=m
+no_defs')

# But the data is in projected space. for a more general script, you would
build this programmatically from the attributes in geo file
wrfproj=CRS("+proj=merc +lat_ts=30 +lon_0=290 +a=6370000.0 +b=6370000.0
+units=m +no_defs")

# create spatial points
# using the extreme corners of the data since raster uses extreme
coordinates. page 3,
https://cran.r-project.org/web/packages/raster/vignettes/Raster.pdf
# there are 16 corner coordinates. 1:4 are cell centers. 13:16 are cell
extremes (see wrf user_guide linked above)
# these coordinates are in longlat so assign geographic system from above
spext <-
SpatialPoints(
list( x=cornerlong[13:16],
y=cornerlat[13:16]),
proj4string = proj1)

# convert to projected space
spext.wrf <- spTransform(spext,wrfproj)

# read the variable wanted
r=raster('../data/geo_em.d01.nc',varname='LANDMASK')
# assign projection detailed in the geo file (defined above as wrfproj)
projection(r)=wrfproj
# change extent to that of the projected spatialpoints
extent(r) <- extent(spext.wrf)
# check the raster.
r
# seem to be missing a couple meters in extent. might be a floating point
issue? you can force the resolution since you know it, but this will change
the bottom right corner of the extent
res(r)=c(dx,dy)
r

        
On Wed, Aug 3, 2016 at 11:51 AM, David Wang <dw2116 at outlook.com> wrote: