netcdf data from Antarctica in rotated latlong - how to reproject to normal geographical coordinate system

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

netcdf data from Antarctica in rotated latlong - how to reproject to normal geographical coordinate system

Matthias Boer-2
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")

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 (direct) | P: +61-(0)2-4570-1941 (admin) | E: [hidden email]

Postal address: Locked Bag 1797 | Penrith, NSW 2751, AUSTRALIA
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: netcdf data from Antarctica in rotated latlong - how to reproject to normal geographical coordinate system

Roger Bivand
Administrator
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

>
> 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 (direct) | P: +61-(0)2-4570-1941 (admin) | E: [hidden email]
>
> Postal address: Locked Bag 1797 | Penrith, NSW 2751, AUSTRALIA
> ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: netcdf data from Antarctica in rotated latlong - how to reproject to normal geographical coordinate system

Michael Sumner-2
On Thu, 5 Apr 2018 at 18:04 Roger Bivand <[hidden email]> 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:
> [hidden email]
> >
> > Postal address: Locked Bag 1797 | Penrith, NSW 2751, AUSTRALIA
> >
> ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
> >
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > 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:
> [hidden email]
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo