[R-sig-Geo] What is the Coordinate Reference System?
Ben Tupper
btupper @ending from bigelow@org
Mon Jun 18 15:22:10 CEST 2018
Hi,
I will blissfully proceed as I have in the past. Thanks to each of you for the clear explanations.
Cheers,
Ben
> On Jun 18, 2018, at 6:26 AM, Michael Sumner <mdsumner using gmail.com> wrote:
>
> Agree with Edzer this application of eqc is anomalous and after it confused
> us initially we now ignore it and apply the unproblematic proj=longlat
> interpretation to these products (after checking).
>
> I might explore the L3 products to find out why this occurred, but the
> ocean colour forum is the usual place to find discussion.
>
> Cheers, Mike
>
> On Mon, 18 Jun 2018, 15:59 Edzer Pebesma, <edzer.pebesma using uni-muenster.de>
> wrote:
>
>> Plate carree is essentially a linear mapping of long and lat to x and y,
>> without changing aspect ratio (scale factor 1 at equator).
>>
>> projecting long/lat towards plate caree using the proj4 string you
>> mention will however change the units from degrees to m. So while the
>> plots of of data in linear long/lat (with aspect ration 1) looks
>> identical to the one obtained in plate caree (with aspect ratio 1), the
>> numerical values along their axes are very different.
>>
>> Many netcdf files carry a grid (raster) with the x and y coordinates of
>> grid cells (centers?) in another coordinate system, but the ones you
>> mention below don't, as far as I can tell. I think what they mean by
>> mentioning plate carree is just saying that this is a regular long/lat
>> degree grid; rather unlucky wording IMO.
>>
>> On 06/15/2018 08:14 PM, Ben Tupper wrote:
>>> 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]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo using r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> Edzer Pebesma
>> Institute for Geoinformatics
>> Heisenbergstrasse 2, 48151 Muenster, Germany
>> Phone: +49 251 8333081
>>
>> _______________________________________________
>> 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]]
More information about the R-sig-Geo
mailing list