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

Thu Jan 5 12:08:28 CET 2017

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.

(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

-- 

Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]