Great example by Michael for this problem
Thanks Michael!
Matt
---------- Forwarded message ----------
From: Michael Sumner <[EMAIL PROTECTED]>
Date: Thu, Jul 31, 2008 at 7:20 PM
Subject: Re: [R-sig-Geo] ncdf file with non-equally spaced grid
To: [EMAIL PROTECTED]
The file really is in Mercator - see the gdalinfo output below - and you
can restore this within a small fudge by reprojecting the bottom corner
and finding an appropriate cell size.
You could use expand.grid(x = lon, y = lat) and work directly with a
SpatialPixelsDataFrame and specify a tolerance for the grid, but that can
be unwieldy for large datasets and I'm more comfortable with the approach
below.
You'll need to check carefully as to whether this is appropriate, but
with the generated TIFF file
you can easily gdalwarp it to longlat if need be.
I think that PROJ.4 will effectively assume your datum is WGS84, but you
should check this with the data source and specify it.
I'm very interested to hear about alternative approaches to this.
HTH
Cheers, Mike.
require(ncdf)
f <- open.ncdf("20071003.276.0237.n17.nc")
lon <- get.var.ncdf(f, "lon")
lat <- get.var.ncdf(f, "lat")
mcsst <- get.var.ncdf(f, "mcsst")
## check the overall alignment
image(lon, lat, mcsst)
## get rgdal for reprojecting and file I/O
library(rgdal)
## get the bottom left corner (not sure if this is cell centre or corner
- assume centre)
cc.offset <- as.vector(project(matrix(c(min(lon), min(lat)), nrow = 1),
"+proj=merc"))
## reproject one column for latitudes
xy.col <- as.matrix(expand.grid(x = min(lon), y = lat))
xy.col.merc <- project(xy.col, "+proj=merc")
## within metres, seems OK
range(diff(xy.col.merc[,2]))
mean(range(diff(xy.col.merc[,2])))
xy.row <- as.matrix(expand.grid(x = lon, y = min(lat)))
xy.row.merc <- project(xy.row, "+proj=merc")
## again, not exact but only within metres
range(diff(xy.row.merc[,1]))
## generate Mercator grid topology
gt <- GridTopology(cc.offset, c(mean(range(diff(xy.row.merc[,1]))),
mean(range(diff(xy.col.merc[,2])))), dim(mcsst))
## flip vertically to get from image() to SGDF alignment
d <- SpatialGridDataFrame(gt, data.frame(mcsst =
as.vector(mcsst[,ncol(mcsst):1])), CRS("+proj=merc"))
## write out to check in GIS (I haven't checked very carefully, but it
seems fine)
writeGDAL(d, "out.tif")
sessionInfo()
R version 2.7.1 (2008-06-23)
i386-pc-mingw32
locale:
LC_COLLATE=English_Australia.1252;LC_CTYPE=English_Australia.1252;LC_MONETARY=English_Australia.1252;LC_NUMERIC=C;LC_TIME=English_Australia.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rgdal_0.5-24 sp_0.9-24 ncdf_1.6
loaded via a namespace (and not attached):
[1] grid_2.7.1 lattice_0.17-8 tools_2.7.1
NETCDF file with subdatasets
C:\temp\DATA\misc\netCDF\R-sig-geo>gdalinfo 20071003.276.0237.n17.nc
Driver: netCDF/Network Common Data Format
Files: 20071003.276.0237.n17.nc
Size is 512, 512
Coordinate System is `'
Metadata:
NC_GLOBAL#Conventions=CF-1.0
NC_GLOBAL#creator_name=Matthew Oliver
[EMAIL PROTECTED]
NC_GLOBAL#institution=Rutgers University Coastal Ocean Observation
Laboratory (RU-COOL) http://rucool.marine.rutgers.edu NC_GLOBAL#url=
http://marine.rutgers.edu/cool/sat_data NC_GLOBAL#source=satellite
radiometer observation NOAA/POES AVHRR instrument
NC_GLOBAL#groundstation=RU-COOL L-band receiver at Rutgers University,
New Jersey, USA
NC_GLOBAL#summary=MCSST calculation and image navigation by TeraScan
software; Regridded to Mercator lon/lat projection for RU-COOL bigbight
subregion
NC_GLOBAL#history=
Subdatasets:
SUBDATASET_1_NAME=NETCDF:"20071003.276.0237.n17.nc":mcsst
SUBDATASET_1_DESC=[1x1222x1183] mcsst (32-bit floating-point)
SUBDATASET_2_NAME=NETCDF:"20071003.276.0237.n17.nc":landmask
SUBDATASET_2_DESC=[1222x1183] landmask (16-bit integer)
SUBDATASET_3_NAME=NETCDF:"20071003.276.0237.n17.nc":cloud_mask
SUBDATASET_3_DESC=[1x1222x1183] cloud_mask (16-bit integer)
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 512.0)
Upper Right ( 512.0, 0.0)
Lower Right ( 512.0, 512.0)
Center ( 256.0, 256.0)
SUBDATASET within the file
C:\temp\DATA\misc\netCDF\R-sig-geo>gdalinfo
NETCDF:"20071003.276.0237.n17.nc":mcsst
Warning 1: Latitude grid not spaced evenly.
Seting projection for grid spacing is within 0.1 degrees threshold.
Driver: netCDF/Network Common Data Format
Files: none associated
Size is 1183, 1222
Coordinate System is `'
Origin = (-77.011923944889588,46.008670050818644)
Pixel Size = (0.011854481576058,-0.009016432203688)
Metadata:
NC_GLOBAL#Conventions=CF-1.0
NC_GLOBAL#creator_name=Matthew Oliver
[EMAIL PROTECTED]
NC_GLOBAL#institution=Rutgers University Coastal Ocean Observation
Laboratory (RU-COOL) http://rucool.marine.rutgers.edu NC_GLOBAL#url=
http://marine.rutgers.edu/cool/sat_data NC_GLOBAL#source=satellite
radiometer observation NOAA/POES AVHRR instrument
NC_GLOBAL#groundstation=RU-COOL L-band receiver at Rutgers University,
New Jersey, USA
NC_GLOBAL#summary=MCSST calculation and image navigation by TeraScan
software; Regridded to Mercator lon/lat projection for RU-COOL bigbight
subregion
NC_GLOBAL#history=
mcsst#units=Celsius
mcsst#missing_value=-9.990000e+002
mcsst#long_name=Multichannel Sea Surface Temperature
mcsst#_FillValue=-9.990000e+002
time#units=seconds since 1970-01-01 00:00:00
time#long_name=Time
lon#units=degrees_east
lon#long_name=Longitude
lat#units=degrees_north
lat#long_name=Latitude
Corner Coordinates:
Upper Left ( -77.0119239, 46.0086701)
Lower Left ( -77.0119239, 34.9905899)
Upper Right ( -62.9880722, 46.0086701)
Lower Right ( -62.9880722, 34.9905899)
Center ( -69.9999981, 40.4996300)
Band 1 Block=1183x1 Type=Float32, ColorInterp=Undefined
NoData Value=-999
Metadata:
NETCDF_VARNAME=mcsst
NETCDF_DIMENSION_time=1191379020
NETCDF_time_units=seconds since 1970-01-01 00:00:00
==============Original message text===============
On Fri, 01 Aug 2008 7:31:49 +1000 "Matt Oliver" wrote:
Sorry, dlat is what changes
an example file can be found at
http://www.ocean.udel.edu/CMS/moliver/20071003.276.0237.n17.nc
use
require(ncdf)
f <- open.ncdf("20071003.276.0237.n17.nc"<
http://www.ocean.udel.edu/CMS/moliver/20071003.276.0237.n17.nc>)
lon <- get.var.ncdf(f, "lon")
lat <- get.var.ncdf(f, "lat")
Thanks in advance
Matt
On Thu, Jul 31, 2008 at 5:03 PM, Michael Sumner <[EMAIL PROTECTED]>wrote:
> You can use gdalwarp in FWTools to re-grid the data with control points
> (perhaps with with some effort), but as you say that will interpolate to a
> regular grid.
>
> Is it definitely dlon that increases? Can you post a link to an example
> file?
>
> Some netCDF files I've seen are actually using Mercator, but store the
> pixel coordinates as lon/lat.
>
> There are totally irregular grids as well though, so as a raster you would
> have to resample it in some way. Otherwise you can interpret the data as a
> SpatialPointsDataFrame and use it as you would in Matlab. Note that the
> image() function will handle irregular (but xy locked) grid spacings, so
you
> can jig a function to plot these things as an image pretty easily.
>
> Cheers, Mike.
>
>
> Matt Oliver wrote:
>
>> Hi All, I have CF compliant ncdf files of sea surface temperature. The
>> lon,
>> lat grid is known but not equally spaced (ie dlon increases). Is there a
>> clean way to get this matrix ported over to a geo-tiff or some other file
>> that Arc will read? I have talked with some arc users that tell me that
>> the
>> grid needs to be equally spaced. I guess I could interpolate to an
equally
>> spaced grid, but I would rather preserve the original grid if possible.
>>
>> Thanks in advance
>>
>> Matt
>>
>>
>>
>>
>
>
--
Matthew J. Oliver
Assistant Professor
College of Marine and Earth Studies
University of Delaware
700 Pilottown Rd.
Lewes, DE, 19958
302-645-4079
http://www.ocean.udel.edu/people/profile.aspx?moliver===========End of
original message text===========
If it wasn't backed-up, then it wasn't important. ~ Anon sysadmin
Mike Sumner (Phd. Candidate)
http://www.antcrc.utas.edu.au/~mdsumner/<http://www.antcrc.utas.edu.au/%7Emdsumner/>
http://www.zoo.utas.edu.au/awru/ IASOS/AWRU
University of Tasmania
Private Bag 80
Hobart Tasmania 7001
AUSTRALIA
Email: [EMAIL PROTECTED]
Phone: 03 6226 1752 (W)
0408599921 (M)
Fax: 03 6226 2745
--
Matthew J. Oliver
Assistant Professor
College of Marine and Earth Studies
University of Delaware
700 Pilottown Rd.
Lewes, DE, 19958
302-645-4079
http://www.ocean.udel.edu/people/profile.aspx?moliver
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
