[R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
Agus Camacho
agus.camacho at gmail.com
Wed Feb 24 00:30:49 CET 2016
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