[R-sig-Geo] reproject coordinates from meters to decimal degrees

Thu Jan 5 12:24:07 CET 2017

On Thu, 5 Jan 2017, Michael Sumner wrote:

> On Thu, 5 Jan 2017 at 19:19 Nikolioudakis Nikolaos <
> nikolaos.nikolioudakis at imr.no> wrote:
>
> Thank you for the prompt response Mike. It would be great if you could do
> so and I can later post the solution on the list in case someone needs it
> in the future
>
>
>
>
> Ugh, it's not what I thought, raster is just wrong on what this file has in
> it, but there's no point laying blame, NetCDF is so general and so commonly
> abused that software everywhere has endless heuristics to figure out what
> to assume.
>
> Basically your data ship with longitude/latitudes, and you need to figure
> out what that is in the native polar stereographic projection as a regular
> grid for use in raster().  You can't easily reproject that raster as it is,
> and I'd advise against it since it's lossy, and doesn't scale well to long
> time series, you could use nn techniques to transform to the longlat values
> in a brute force way, but that's a horrible hack - I think it's worth
> figuring this one out.
>
> I've put my startings here: https://github.com/mdsumner/ranch/issues/3
>
> If I can I'll figure it out and let you know.  Feel free to ping on that
> issue, but also share more info you find from the sources -  it's a
> repository I made to collect NetCDF adventures - it  might be worthwhile .
> . .
>
> Similar data (from raw binary but using the same-mostly projection) is
> wrapped by us here:
> https://github.com/AustralianAntarcticDivision/raadtools/blob/master/R/ice.R
> This is also built into raster itself (but only to read those particular
> binary files, not NetCDF ones).
>
> I kind of agree with Roger about a vignette for rgdal, but really the
> information belongs at the source of the data (and inside the NetCDF file
> itself), and I'd be more interested if rgdal was easier to work with as a
> collaborator. (Most users would think it belongs in raster, since that
> hides rgdal for the most part - neither package is easy to contribute to).
>
> I agree we need a better community resource for these file-bogeys, but
> it's something we share with the GDAL community and many others so it's
> bigger than any particular software user group.

Yes, I agree. Does the use of +proj=obtran - oblique transformations - 
come into the same general area of file-bogeys which could be documented 
in the same place?

Roger

>
> (For this problem generally, we just collect our fixes in our local tool
> kit (raadtools) for files that will be of long-term use, and because it
> provides the rest of the workflows we need on top of huge file collections
> on local shares. I don't see any general solutions that are better than
> that on the horizon, though incremental improvements to raster, GDAL and
> NetCDF have seen a lot of our workarounds updated).
>
> Cheers, Mike.
>
>
>
> Cheers,
>
> Nikos
>
>
>
> *From:* Michael Sumner [mailto:mdsumner at gmail.com]
> *Sent:* Thursday, January 05, 2017 9:01 AM
> *To:* Nikolioudakis Nikolaos <nikolaos.nikolioudakis at imr.no>;
> r-sig-geo at r-project.org
> *Subject:* Re: [R-sig-Geo] reproject coordinates from meters to decimal
> degrees
>
>
>
> This data is in polar stereographic so it's best to ignore the lonlats and
> apply the right affine transform (extent and resolution) directly. I.e. no
> transformation or hacky workarounds required. We do this routinely from
> similar files, happy to explain off-list.
>
> For annoying reasons this is poorly understood on both the GIS and NetCDF
> (modelling) sides.
>
> Cheers, Mike
>
>
>
> On Thu, 5 Jan 2017, 18:43 Nikolioudakis Nikolaos <
> nikolaos.nikolioudakis at imr.no> wrote:
>
> Dear members of the list,
> My question might be trivial but please bear with me as I am quite new in
> GIS and R so my terminology and/or expressions might not be the right ones.
> I have a set of coordinates in lat/lon for which I need to extract values
> from a netCDF file. The file is downloaded from http://marine.copernicus.eu/
> (needs registration to reach data from there). For this reason you can find
> in the following link a .nc subset (single raster) from a multiband .nc I
> am using.
> .nc file:
>
> https://dl.dropboxusercontent.com/u/2156982/subset.nc
>
> The coordinates of my locations (in lon/lat) are in the coords.csv file can
> be found here:
>
> https://dl.dropboxusercontent.com/u/2156982/coords.csv
>
> The raster's projection is (according to the data source EPSG: 3411 (NSIDC
> Polar Stereographic North)) so the obvious thing to do is reproject my
> points to the same projection as the raster's in order to overplot them
> (and subsequently extract values from the raster based on the points).
> Reprojecting the rasters would also be an option, but in my specific case I
> have multiband (many dates) netCDFs and also many different variables in
> each netCDF so reprojecting the coordinates if possible and performing the
> extraction from the epsg 3411 rasters would be ideal for me. My problem is
> that when I reproject my coordinates I get them in meters and not in
> decimal degrees as the raster is.
>
>
>
> So the question is how can I obtain my reprojected coordinates in decimal
> degrees instead of meters in order to be able to perform the extraction of
> values from the raster (I have done this for another dataset that both the
> locations and the rasters were in the same projection and units). I
> understand that my problem is probably related with how to declare that the
> output units of my reprojected coordinates should be in dd instead of
> meters but I am not sure on the syntax of the proj4 string and on the web
> the inverse procedure is usually described (dd to meters). The r code I am
> using is the following:
>
> library(raster)
>
> library(rgdal)
>
>
>
> #setwd to the location of the files
>
> subset.r = raster('subset.nc', varname='chla')
>
> #print(subset.r)
>
>
>
> locs <- read.table('coords.csv',header = T,sep = ';')
>
> coordinates(locs) <- ~ LON+LAT
>
> # par(mfrow=c(1,2))
>
> # plot(subset.r)
>
> # plot(locs at coords)
>
>
>
> #set the projection to EPSG: 4326 for the locs
>
> wgs84 <- CRS('+init=epsg:4326')
>
> proj4string(locs) <- wgs84
>
> #reproject locs to polar stereographic (EPSG 3411)
>
> locs.repr <- spTransform(locs, crs(subset.r))
>
> #compareCRS(locs.repr, subset.r)
>
> plot(subset.r)
>
> plot(locs.repr at coords<mailto:locs.repr at coords>)
>
> # extent(locs.repr)
>
> # extent(subset.r)
>
>
>
> Any help would be greatly appreciated as well as any comments on things I
> might be forgetting to provide so as to make my problem better understood.
>
> Best regards,
> Nikos
> __________________________
> Nikos Nikolioudakis
> PostDoc Researcher
> Pelagic Fish Research Group
> Institute of Marine Research (IMR)
> P.O. Box 1870 - Nordnes
> 5817 Bergen, Norway
> E-mail:nikolaos.nikolioudakis at imr.no
> Tel:0047489981 <0047489981>
>
> [cid:image001.jpg at 01D0B00F.4BA30FC0]
> __________________________
>
>        [[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; 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