Agus,
sorry for the addon but I think I have to provide a correction of the
corner coordinates (e.g. the extent values):
In the example that you have posted below I did calculate the extent
using the domain center coordinate and the WRF grid resolution in meter
and the number of rows and cols.
Since Dominik provides a link to the file description of the netcdf file
I think it is more accurate to reproject the corner coordinates as
given by the netcdf header variables (NC_GLOBAL#corner_lats,
NC_GLOBAL#corner_lons). Assuming that your variable "dctmx" (which I can
not identfy it in the nc file) is of type "Mass" (staggered = M) the
correct corner coordinates are stored as the first 4 entrys of the dump
snippet below:
lats<-"
NC_GLOBAL#corner_lats={12.355667,50.26619,50.26619,12.355667,12.308136,50.18787,50.18787,12.308136,12.210785,50.403816,50.403816,12.210785,12.163345,50.325382,50.325382,12.163345}"
lons<-"
NC_GLOBAL#corner_lons={-131.43678,-151.29639,-48.703613,-68.563232,-131.5851,-151.51157,-48.488434,-68.414917,-131.38828,-151.41891,-48.581085,-68.611725,-131.53641,-151.63441,-48.36557,-68.463593}"
after cleaning and converting the strings you may calculate the corner
coordinates:
library(proj4)
## project mass corner coordinates to lambertian
llMass <- ptransform(cbind(clon[1],clat[1])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
ulMass <- ptransform(cbind(clon[2],clat[2])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
lrMass <- ptransform(cbind(clon[3],clat[3])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
urMass <- ptransform(cbind(clon[4],clat[4])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
wrfLccExtMass<-extent(ulMass[1],lrMass[1],llMass[2],ulMass[2])
According to this the correct extent for mass variables should be:
extent(wrfLccExtMass)
class : Extent
xmin: -3575343
xmax: 3575342
ymin: -2293058
ymax: 2306330
hope this is correct now
cheers Chris
Am 24.02.2016 um 05:12 schrieb Agus Camacho:
With the help of everybody i finally got this to work, so here is a script
that does the job of reprojecting both, a raster layer obtained from a .nc
and some locations in order to overplot them using plot.raster or mapview.
Im using a combination of the advices of Dominik, Michael and Chris.
require(ncdf4)
require(raster)
setwd("C:/~")
r=raster("geo_em.d01.nc",
varname="dctmx")# days of ctmax events
# Set extent and projections of rasters for plotting
# chris gave me the orig data fom the nc file because i could not install
gdal
xmin= -3545999
xmax= 3546000
ymin = -2286000
ymax=2286000
pr<- "+proj=lcc +lat_1=25 +lat_2=45 +lat_0=38.08 +lon_0=-100 +x_0=0
+y_0=0"
wrfLccExt<-extent(xmin,xmax,ymin,ymax)
extent(r) <- extent(wrfLccExt)
projection(r) <- pr
# get and prepare urosaurus locations
x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
x=x[,1:3]
colnames(x)
coordinates(x)=~lon+lat
proj4string(x)<- CRS("+proj=longlat +datum=WGS84")
x=spTransform(x, pr)
# get prepare and plot wrld_simpl
require(maptools)
require(sp)
data(wrld_simpl)
namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
"Canada", "Mexico")), prj)
#plot with raster
plot(r, cex.main=.7,legend=F)
points(x)
plot(namer, add = TRUE)
# plot with mapview (cool!)
m=mapview(r)
u=mapview(x)
m+u
Thanks to all!
Agus
2016-02-23 13:11 GMT-07:00 Agus Camacho :
Thanks for that Dominik,
Giving that projection to either the locations, the raster layer generated
from the .nc file, or both, still did not work. I keep having locations
that should be on land falling far on the sea. Might this be a problem
derived from using raster with a file whose original grid distances are not
constant?
Here is a link with the original file which has the original coordinate
data.
https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
2016-02-23 12:07 GMT-07:00 Dominik Schneider <
dominik.schnei...@colorado.edu>:
This looks like WRF data. I just
dealt with this.
The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
+a=637 +b=637 +units=m
+proj=lcc which is usually what wrf is run with.
The tricky part is:
+lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
because every WRF run is different (the WRF Preprocessing System
optimizes the projection for the domain).
and then there is probably no shift so you need(?) +x_0=0 +y_0=0
This gives:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
+ellps=sphere +a=637 +b=637 +units=m +no_defs
But, wrf users like to give out lat and long so you need to assign it:
+proj=longlat +ellps=sphere +a=637 +b=637 +units=m +no_defs
and then reproject the lat/long to lcc coordinates using this string:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
+ellps=sphere +a=637 +b=637 +units=m +no_defs
One word of c
