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

Michael Sumner mdsumner at gmail.com
Wed Feb 24 00:24:46 CET 2016 

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

	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list