[R-sig-Geo] Set up raster data in WRF Mercator projection
Dominik Schneider
dosc3612 at colorado.edu
Wed Aug 3 22:54:35 CEST 2016
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:
> 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]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list