[R-sig-Geo] netcdf data from Antarctica in rotated latlong - how to reproject to normal geographical coordinate system

[Michael Sumner]

On Thu, 5 Apr 2018 at 18:04 Roger Bivand <Roger.Bivand at nhh.no> wrote:

> On Thu, 5 Apr 2018, Matthias Boer wrote:
>
> > Hi all
> >
> > I am having trouble extracting some climate time series for a series of
> field sites in Antarctica from a netcdf file with a rotated pole
> projection. The file (gridded monthly temperatures for 1979-2016) includes
> these details about the projection:
> > float rotated_pole[]
> >            grid_mapping_name: rotated_latitude_longitude
> >            grid_north_pole_latitude: -180
> >            grid_north_pole_longitude: -150
>
> >            proj4_params: -m 57.295779506 +proj=ob_tran +o_proj=latlon
> > +o_lat_p=-180.0 +lon_0=30.0
>
> >            proj_parameters: -m 57.295779506 +proj=ob_tran +o_proj=latlon
> > +o_lat_p=-180.0 +lon_0=30.0
>
> >            projection_name: rotated_latitude_longitude
> >            long_name: projection details
> >            EPSG_code:
> >
> > I can read the netcdf file into R:
> > tst.s <- stack("filename.nc", varname="var1")
>
> You'll need to look for the arguments supporting +proj=ob_tran
> (use_ob_tran=) in rgdal::spTransform(). +proj=ob_tran is a complete mess,
> and until the next release of GDAL, the GDAL netcdf driver will not
> understand these either. Your best bet may be not to warp the input grids,
> but to transform the field sites point coordinates to ob_tran to extract
> the data you need (checking that the transformation seems sensible).
>
> Roger
>
>
Just to echo Roger's expert knowledge of the vagaries of ob_tran support
here - this is a fraught area, and the general advice to *transform your
query to the grid, rather than vice versa* is so very important. If you can
make a map with the grid - plot(tst.s[[1]]) - and then add vector (points,
lines, polygons) data to that  by transforming them to that space then you
know for sure it can all work, and be optimized and automated.The effort is
in getting the transformation of a few points assured, and then you can
bundle it all up.

It's not obvious, and I'm happy to discuss in detail offline to get a
solution - it always takes a bit of back and forth and experimenting -
which is hard to get across in a public forum - for my own purposes, seeing
real world example data (especially Antarctic!) is always welcome as grist
for understanding.

Can you share a file? (Offline if need be).

Cheers, Mike.


>
> > but this gives a couple of warnings that raster can't understand the
> projection:
> > Warning messages:
> > 1: In .stackCDF(x, varname = varname, bands = bands) :
> >  tskin has 4 dimensions, I am using the last one
> > 2: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
> > grid_north_pole_latitude=-180; grid_north_pole_longitude=-150;
> proj4_params=-m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lat_p=-180.0
> +lon_0=30.0; proj_parameters=-m 57.295779506 +proj=ob_tran +o_proj=latlon
> +o_lat_p=-180.0 +lon_0=30.0; projection_name=rotated_latitude_longitude;
> long_name=projection details; EPSG_code=
> > 3: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
> > grid_north_pole_latitude=-180; grid_north_pole_longitude=-150;
> proj4_params=-m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lat_p=-180.0
> +lon_0=30.0; proj_parameters=-m 57.295779506 +proj=ob_tran +o_proj=latlon
> +o_lat_p=-180.0 +lon_0=30.0; projection_name=rotated_latitude_longitude;
> long_name=projection details; EPSG_code=
> >
> >
> > Could anyone point me to some code to assign the correct CRS to the
> > imported stack 'tst.s' and then to project it to the 'normal'
> > geographical coordinate system, so I can extract time series for the
> > field site coordinates which are in 'normal' latlong.
> >
> > Thanks a lot for your time.
> >
> > Cheers,
> > Matthias
> >
> >
> >
> >
> ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
> > Dr. Matthias M Boer
> > Hawkesbury Institute for the Environment<
> http://www.westernsydney.edu.au/hie> | Western Sydney University |
> Hawkesbury Campus, Richmond, NSW 2753, AUSTRALIA
> > P: +61-(0)2-4570-1373 <+61%202%204570%201373> (direct) | P:
> +61-(0)2-4570-1941 <+61%202%204570%201941> (admin) | E:
> m.boer at westernsydney.edu.au
> >
> > Postal address: Locked Bag 1797 | Penrith, NSW 2751, 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 <+47%2055%2095%2093%2055>; e-mail:
> Roger.Bivand at nhh.no
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
> _______________________________________________
> 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]]