[R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Dominik Schneider Dominik.Schneider at colorado.edu
Wed Feb 24 05:09:18 CET 2016 

just glancing through the email chain. looks like you were able to get
things plotted....but you are still having an issue with some points?

here is a link for information regarding the grid staggering used in WRF.
http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm

> Some remarks:
> (1) I just took the first pair of  coordinates as derived from gdalinfo("
> geo_em.d01.nc")
> you will find 4 different coordinate pairs (i did not proof which one is
> right
> The data is staggered (as outlined by Dominik) So some of the corner
> coordinates belongs to the staggered data and the others coordinates  to
> the unstaggered ones.



On Tue, Feb 23, 2016 at 4:30 PM, Agus Camacho <agus.camacho at gmail.com>
wrote:

> Thanks Michael, that was a good job and despite the map looks nice now,
> the locations cant still be plotted adequately, either as raw lon lat or as
> spatial points with different projections.
>
>
> 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
> projection(x)="+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
> x=spTransform(x, CRS("+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"))
> points(x)
>
> 2016-02-23 16:24 GMT-07:00 Michael Sumner <mdsumner at gmail.com>:
>
>>
>>
>> On Wed, 24 Feb 2016 at 06:07 Dominik Schneider <
>> Dominik.Schneider at 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.camacho at 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 <mdsumner at gmail.com>:
>>>>
>>>> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <Roger.Bivand at 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
>>>> > > >> frtog at 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-bounces at r-project.org] On
>>>> Behalf Of
>>>> > > Agus Camacho
>>>> > > >> Sent: 22. februar 2016 19:20
>>>> > > >> To: tech at 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_dev at 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 at 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.Bivand at 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 at 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 at 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 at 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
>>
>>
>
>
> --
> 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]]

More information about the R-sig-Geo mailing list