[R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
Chris Reudenbach
reudenbach at uni-marburg.de
Wed Feb 24 09:36:22 CET 2016
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.
>>
>
>
More information about the R-sig-Geo
mailing list