adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
On Tue, 23 Feb 2016, Agus Camacho 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
Thanks for the link. The proj4 string is missing +units=km (the spherical axes are certainly not in metres), but is wrong anyway. This says:
gI <- GDALinfo("results_us_future_output_none_0.nc")
Error in .local(.Object, ...) :
No UNIDATA NC_GLOBAL:Conventions attribute
The values of the coordinate variables are:
library(ncdf4)
r <- nc_open("results_us_future_output_none_0.nc")
pts <- cbind(c(ncvar_get(r, "lon")), c(ncvar_get(r, "lat")))
summary(pts)
V1 V2 Min. :-151.30 Min. :12.36 1st Qu.:-120.44 1st Qu.:25.45 Median :-100.00 Median :35.72 Mean :-100.00 Mean :35.71 3rd Qu.: -79.56 3rd Qu.:45.91 Max. : -48.70 Max. :58.33 which are obviously not a grid: pts1 <- SpatialPoints(pts) plot(pts1) # gridded(pts1) <- TRUE fails miserably This is definitely not lcc, irrespective of what the originator says. Mike - any forensics? Roger
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:
and a wrld_simpl object: library(maptools) data(wrld_simpl) Over this raster layer
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
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