[R-sig-Geo] Set up raster data in WRF Mercator projection

[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

________________________________
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<mailto: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<http://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<http://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]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:R-sig-Geo at r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


	[[alternative HTML version deleted]]