On Wed, 24 Feb 2016 at 06:07 Dominik Schneider < dominik.schnei...@colorado.edu> wrote:
> This looks like WRF <http://www.wrf-model.org/index.php> data. I just > dealt with this. > The data is on a sphere as opposed to WGS84 so you need +ellps=sphere > +a=6370000 +b=6370000 +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=6370000 +b=6370000 +units=m +no_defs > > This looks right to me: library(raster) library(rgdal) prj <- "+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs" r <- raster("results_us_future_output_none_0.nc", varname = "dctmx") ## use print(r) to see the nc header dump ## lons/lats from NetCDF lon <- raster("results_us_future_output_none_0.nc", varname = "lon") lat <- raster("results_us_future_output_none_0.nc", varname = "lat") ## the true projected coordinates (??) xy <- project(cbind(values(lon), values(lat)), prj) ## looks right plot(xy, pch = ".") ## this fails, so it's not exactly right ##grd <- rasterFromXYZ(cbind(xy, 0) ## could use sp::points2grid tools to control tolerances, or just fudge it ex <- extent(xy) resolution <- c(round(max(diff(sort(unique(xy[,1]))))), round(max(diff(sort(unique(xy[,2])))))) ex <- ex + c(-1, 1, -1, 1) * rep(resolution / 2, 2) ## finally mp <- setExtent(r, ex) projection(mp) <- prj library(maptools) data(wrld_simpl) namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States", "Canada", "Mexico")), prj) plot(mp) plot(namer, add = TRUE) ("Why does this crucial metadata get thrown away?", he sighs . . .) HTH > But, wrf users like to give out lat and long so you need to assign it: > +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +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=6370000 +b=6370000 +units=m +no_defs > > One word of caution, make sure you received the correct coordinates. Some > variables are run cell center while some are run at cell edge. It looks > like from your .nc file it was made by your collaborator so I assume they > are right. > > That said, another word of caution, I found that the XLAT and XLONG > variables from WRF output aren't very precise. There is a "geogrid" file > from the preprocessing system that has the domain corners, resolution, nrow > and ncol from which you can make a better grid using the native projection > system (in my case it was a 4km grid). I suggest you try to get those. > > I hope this helps... I have to run but wanted to save people too much head > scratching. I can get you running with more help tonight if you need. > Dominik > > > On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <agus.cama...@gmail.com> > wrote: > >> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring >> more giant crab this time. >> >> The creator of the .nc file also looked at this webpage: >> http://www.pkrc.net/wrf-lambert.html >> It seemed like the right proj4 string might be this one: >> >> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 >> +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs >> >> However this projection also does not allow me to adequately plot the >> locations on the raster. >> >> Here is the .nc file. it contains several layers. >> >> >> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0 >> >> >> >> 2016-02-23 2:25 GMT-07:00 Michael Sumner <mdsum...@gmail.com>: >> >> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <roger.biv...@nhh.no> wrote: >> > >> > > On Tue, 23 Feb 2016, Alex Mandel wrote: >> > > >> > > > I made an attempt at it too. Investigating the original data, I'm >> not >> > > > sure that the projection information supplied is correct for the >> data >> > > > linked. When I load up the data in a unprojected space, the >> coordinates >> > > > don't look at all similar to any Lambert projected data I have, they >> > > > actually look like Lat/Lon in some unprojected coordinate system, >> > > > perhaps a different spheroid than expected. >> > > >> > > Does anyone have a link to the original data? Is is possible that >> this is >> > > the General Oblique Transformation used by modellers - that is >> something >> > > that feels like longlat but is recentred and oblique? Example at the >> very >> > > end of my GEOSTAT talk last year (slides 81-83): >> > > >> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf >> > > >> > > Roger >> > > >> > > >> > For what it is worth, the General Oblique Transformation is not the only >> > example - it's very common for modellers to have a mesh that has the >> > "mostly-properties" of a projection, but is not actually describable >> with >> > standard transform + affine parameters. The main cases that I've seen >> are >> > polar stereographic, equal area or oblique Mercator. Often they really >> are >> > simple transforms and you can reconstruct without loss, but it's not >> > usually possible to tell without exploration. It's an interesting >> > dis-connect to see code that builds a mesh with certain properties, then >> > only stores longitudes and latitudes - when it could be done with >> standard >> > tools and be stored and used much more efficiently. >> > >> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area >> > terminology conflated in this context too. ) >> > >> > I'm also interested to explore the original data. >> > >> > Cheers, Mike. >> > >> > >> > >> > > > >> > > > -Alex >> > > > >> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote: >> > > >> Hi >> > > >> >> > > >> I tried to make it work but I had to give up. I wanted to reproject >> > the >> > > Lamberth conformal conic coordinates to long-lat but it didn't work. >> > > >> >> > > >> Perhaps someone can see what I did wrong. Here is what I did (data >> in >> > R >> > > binary format and figure in png format both attached): >> > > >> >> > > >> library(raster) >> > > >> library(maptools) >> > > >> data(wrld_simpl) >> > > >> >> > > >> r <- raster("raster.grd") >> > > >> projection(r) >> > > >> ## > NA >> > > >> >> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = >> ",") >> > > >> coordinates(uro) <- ~lon+lat >> > > >> >> > > >> ## Set projections for the 3 data sets >> > > >> >> > > >> ## Lamberth's confocal conic projection with given parameters >> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 >> +lat_2=45.0 >> > > +ellps=WGS84" >> > > >> projection(r) >> > > >> >> > > >> ## Assume that lon, lat are geographical coordinates (degrees >> decimal) >> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84") >> > > >> >> > > >> ## wrld_simpl is in geographical coordinates >> > > >> proj4string(wrld_simpl) >> > > >> >> > > >> ## Make figure in png format >> > > >> ## Of course plotting data with 2 different projections will give >> > > >> ## some distortions >> > > >> pdf("uro.png") >> > > >> >> > > >> plot(r) >> > > >> points(uro) >> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of >> 'r' >> > > >> >> > > >> dev.off() >> > > >> >> > > >> extent(r) >> > > >> ## class : Extent >> > > >> ## xmin : -131.4368 >> > > >> ## xmax : -68.56323 >> > > >> ## ymin : 12.35567 >> > > >> ## ymax : 50.26619 >> > > >> >> > > >> ## Reproject the raster to long-lat >> > > >> ## This doesn't work (collapsed domain) >> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs >> > > +ellps=WGS84 +towgs84=0,0,0") >> > > >> >> > > >> ## Because >> > > >>> extent(rp) >> > > >> ## class : Extent >> > > >> ## xmin : -100.0015 >> > > >> ## xmax : -99.68557 >> > > >> ## ymin : 37.70658 >> > > >> ## ymax : 38.00046 >> > > >> >> > > >> ## Save data in R binary format >> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData") >> > > >> >> > > >> >> > > >> Yours sincerely / Med venlig hilsen >> > > >> >> > > >> Frede Aakmann Tøgersen >> > > >> Specialist, M.Sc., Ph.D. >> > > >> Plant Performance & Modeling >> > > >> >> > > >> Technology & Service Solutions >> > > >> T +45 9730 5135 >> > > >> M +45 2547 6050 >> > > >> fr...@vestas.com >> > > >> http://www.vestas.com >> > > >> >> > > >> Company reg. name: Vestas Wind Systems A/S >> > > >> This e-mail is subject to our e-mail disclaimer statement. >> > > >> Please refer to www.vestas.com/legal/notice >> > > >> If you have received this e-mail in error please contact the >> sender. >> > > >> >> > > >> >> > > >> >> > > >> -----Original Message----- >> > > >> From: R-sig-Geo [mailto:r-sig-geo-boun...@r-project.org] On >> Behalf Of >> > > Agus Camacho >> > > >> Sent: 22. februar 2016 19:20 >> > > >> To: t...@wildintellect.com >> > > >> Cc: r-sig-geo >> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a >> > > reference system implicit in a .nc file >> > > >> >> > > >> Thanks Alex, but the locations still fall in the sea when i plot >> them >> > > using >> > > >> your recommended Solution. I looked at the sites you proposed and >> they >> > > have >> > > >> other values for lat_1, lat_0, etc.. >> > > >> >> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <tech_...@wildintellect.com>: >> > > >> >> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote: >> > > >>>> Dear all, >> > > >>>> >> > > >>>> Im trying to overlap these points: >> > > >>>> >> > > >>> >> > > >> > >> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0 >> > > >>>> >> > > >>>> and a wrld_simpl object: >> > > >>>> library(maptools) >> > > >>>> data(wrld_simpl) >> > > >>>> >> > > >>>> Over this raster layer >> > > >>>> >> > > >>> >> > > >> > >> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0 >> > > >>>> >> > > >>>> This rastr comes from a .nc file without a reference system. The >> > > author >> > > >>> of >> > > >>>> that .nc file gave me the following data about the .nc. >> > > >>>> >> > > >>>> The projection is *Lambert conformal conic* projection >> > > >>>> CEN_LAT = 38.0 >> > > >>>> CEN_LON = -100.0 >> > > >>>> TRUELAT1 = 25. >> > > >>>> TRUELAT2 = 45. >> > > >>>> >> > > >>>> However, despite i have gone through many sites in the internet, >> i >> > > cant >> > > >>>> figure it out: >> > > >>>> >> > > >>>> a) if that is all the data i need to set a reference system for >> my >> > > points >> > > >>>> and the wrld_simp object. >> > > >>>> >> > > >>>> b) how to change a typical CRS object with such data >> > > >>>> >> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84") >> > > >>>> >> > > >>>> Where do i enter the TRUELAT and CENLAT values? >> > > >>>> Are there any site that explains easily what the fields in the >> CRS >> > > mean >> > > >>> and >> > > >>>> how to change them? >> > > >>>> >> > > >>>> Thanks in advance. >> > > >>>> >> > > >>> >> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms >> > > >>> https://trac.osgeo.org/proj/wiki/GenParms >> > > >>> >> > > >>> I believe: >> > > >>> +lat_0 = CEN_LAT Latitude of origin >> > > >>> +lat_1 = TRUELAT1 Latitude of first standard parallel >> > > >>> +lat_2 = TRUELAT2 Latitude of second standard parallel >> > > >>> +lon_0 = CEN_LON Central meridian >> > > >>> >> > > >>> proj strings are defined by the proj4 libary. It's website listed >> > above >> > > >>> and the associated mailing lists or gis stackexchange would be the >> > > >>> places to get help on it. >> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs >> > > >>> >> > > >>> It often helps to browse similar projections on >> > > >>> http://spatialreference.org/ >> > > >>> http://epsg.io/ >> > > >>> >> > > >>> Enjoy, >> > > >>> Alex >> > > >>> >> > > >>> >> > > >> >> > > >> >> > > > >> > > > _______________________________________________ >> > > > R-sig-Geo mailing list >> > > > R-sig-Geo@r-project.org >> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > >> > > -- >> > > Roger Bivand >> > > Department of Economics, Norwegian School of Economics, >> > > Helleveien 30, N-5045 Bergen, Norway. >> > > voice: +47 55 95 93 55; fax +47 55 95 91 00 >> > > e-mail: roger.biv...@nhh.no >> > > http://orcid.org/0000-0003-2392-6140 >> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en >> > > http://depsy.org/person/434412 >> > > _______________________________________________ >> > > R-sig-Geo mailing list >> > > R-sig-Geo@r-project.org >> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > >> > -- >> > Dr. Michael Sumner >> > Software and Database Engineer >> > Australian Antarctic Division >> > 203 Channel Highway >> > Kingston Tasmania 7050 Australia >> > >> > [[alternative HTML version deleted]] >> > >> > _______________________________________________ >> > R-sig-Geo mailing list >> > R-sig-Geo@r-project.org >> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > >> >> >> >> -- >> Agustín Camacho Guerrero. >> Doutor em Zoologia. >> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de >> Biociências, USP. >> Rua do Matão, trav. 14, nº 321, Cidade Universitária, >> São Paulo - SP, CEP: 05508-090, Brasil. >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > -- Dr. Michael Sumner Software and Database Engineer Australian Antarctic Division 203 Channel Highway Kingston Tasmania 7050 Australia [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo