Set up raster data in WRF Mercator projection

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

------------------------------
*From:* Dominik Schneider <dosc3612 at colorado.edu>
*Sent:* Wednesday, August 3, 2016 12:52 PM
*To:* David Wang
*Cc:* r-sig-geo
*Subject:* Re: [R-sig-Geo] Set up raster data in WRF Mercator projection

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

From my notes, the mercator WRF projection would look like this:
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:

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

        [[alternative HTML version deleted]]