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

Jonathan Greenberg jgrn at illinois.edu
Fri Feb 26 16:17:04 CET 2016


Folks:

Just noticed this thread -- I see I didn't include a "--config" option with
any of the gdalUtils functions (it isn't one of the documented parameters
on the individual utility website, but it seems it would have allowed you
to run the GDAL_NETCDF_BOTTOMUP without setting a system environment
variable -- see
http://lists.osgeo.org/pipermail/gdal-dev/2014-July/039452.html).  If I
make this tweak and pushed it to r-forge, would one of you be willing to
see if it solves the problem?

I assume this will be something that would be needed for any gdal utility
that allows an "of" to be set?

--j

On Wed, Feb 24, 2016 at 2:37 AM Chris Reudenbach <reudenbach at uni-marburg.de>
wrote:

> 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.000008 +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 <agus.camacho at gmail.com>:
> >
> >> 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.Schneider at colorado.edu>:
> >>
> >>> 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
> >>>
> >>> 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
> >>>>
> >>>
> >>
> >> --
> >> 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.
> >>
> >
> >
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list