[R-sig-Geo] What is the Coordinate Reference System?

[Ben Tupper]

Hi,

That's an interesting approach.  

On the other hand, sometimes I get bitten (well, muddled really) when I use raster to determine the crs of a ncdf resource.  It all sort of depends.

Below is an example for Aqua MODIS imagery from OBPG (https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/Mapped/Daily/9km/sst/2017/ <https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/Mapped/Daily/9km/sst/2017/>).  The OBPG docs (https://oceancolor.gsfc.nasa.gov/products/ <https://oceancolor.gsfc.nasa.gov/products/>) state that the data is served up as 'Plate Carrée' which I think has a crs ...

'+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ' 

... as described here http://spatialreference.org/ref/epsg/32662/ <http://spatialreference.org/ref/epsg/32662/>

You can see below that raster comes up with...

'+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0' 

... which is http://spatialreference.org/ref/epsg/4326/ <http://spatialreference.org/ref/epsg/4326/>

Explicitly setting the crs while reading doesn't help.  I have never resolved if it is an OBPG documentation bug, my confusion about projections, or some other mystery.  So, my usual practice is to assume that raster makes the right call.

#### start

library(raster)
URI = 'https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A2017175.L3m_DAY_SST_sst_9km.nc'

# ~3.8MB
ok = download.file(URI, basename(URI))

R = raster::raster(basename(URI), varname = 'sst')
R
# class       : RasterLayer
# dimensions  : 2160, 4320, 9331200  (nrow, ncol, ncell)
# resolution  : 0.08333334, 0.08333334  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# data source : /Users/Shared/code/R/others/record/spnc/A2017175.L3m_DAY_SST_sst_9km.nc
# names       : Sea.Surface.Temperature
# zvar        : sst

S = raster::raster(basename(URI), varname = "sst", crs = '+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
# Warning message:
# In .getProj(prj, crs) :
#  argument "crs" ignored because the file provides a crs

#### end


Cheers,
Ben


> On Jun 14, 2018, at 6:42 PM, Michael Sumner <mdsumner using gmail.com> wrote:
> 
> Try raster::projection(raster::raster(a)) first, then delve if necessary
> into the ncdump summary with print(raster::raster(a)))
> 
> But, raster(a) may be enough here? To get that stuff more directly from the
> NetCDF API you have to use the inq and att functions of ncdf4 or RNetCDF
> and that's no fun.
> 
> I'd also try raster(rgdal::readGDAL(a)) which is a totally different but
> equivalent read path.
> 
> Cheers, Mike
> 
> On Thu, 14 Jun 2018, 11:03 Sergio Ibarra, <zergioibarra using gmail.com> wrote:
> 
>> Greetings,
>> 
>> I'm want to read NetCDF spatial data from a meteorological model (WRF
>> model). The NetCDF can have one of these three projections:
>> 
>> Polar stereographic
>> Lambert conformal
>> Mercator projection
>> 
>> (See figure
>> 
>> https://user-images.githubusercontent.com/27447280/41389756-ed23f91a-6f68-11e8-961e-6ba901913c54.png
>> )
>> 
>> What is the  Coordinate Reference System for each case?
>> 
>> Thanks!
>> 
>> Here an example with ncdf4 and raster:
>> 
>> library(ncdf4)
>> library(raster)
>> a <- tempfile()
>> download.file(url = "
>> http://www.unidata.ucar.edu/software/netcdf/examples/wrfout_v2_Lambert.nc
>> ",
>>              destfile = a) #78 Mb
>> wrf <- nc_open(a)
>> paste("The file has",wrf$nvars,"variables")
>> paste("The file has",names(wrf$var))
>> lat <- ncvar_get( wrf, "XLAT" )
>> lon <- ncvar_get( wrf, "XLONG" )
>> t2 <- ncvar_get( wrf, "T2" ) - 273.15
>> dim(t2) #73 60 13
>> #image
>> image(t2[, , 12])
>> # raster
>> rt2 <- raster(t(t2[1:dim(t2)[1],
>>                   dim(t2)[2]:1,
>>                   12]),
>>              xmn = min(lon),
>>              xmx = max(lon),
>>              ymn = min(lat),
>>              ymx = max(lat),
>> *              crs="+init=epsg:4326") # ???*
>> spplot(rt2)
>> 
>> 
>> 
>> --
>> ​Sergio Ibarra Espinosa
>> Post Doc
>> PhD in Atmospheric Sciences
>> Instituto de Astronomia, Geofísica e Ciências Atmosféricas
>> Universidade de São Paulo
>> Rua do Matão, 1226
>> São Paulo-SP - Brasil -
>> 05508-090
>> +55-11-97425-3791
>> Skype: sergio_ibarra1
>> 
>>        [[alternative HTML version deleted]]
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo using 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 using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

Ecological Forecasting: https://eco.bigelow.org/






	[[alternative HTML version deleted]]