Re: [R-sig-Geo] ncdf file with non-equally spaced grid

Wed, 06 Aug 2008 02:56:45 -0700 

With FWTools installed, you can do this with the resulting out.tif

gdalwarp out.tif ll.tif -t_srs "+proj=longlat"
That will reproject the result to longlat. I suggest you check the result carefully, and do find out the datum you should use.
I've realized that trying to extract the stuff from the netCDF directly as I do below results in the image being flipped vertically, because of the way it is stored internally. I don't know how to do this otherwise without using the VRT driver.
I'll try to put together an example, but I suggest you explore the command line utilities for GDAL - most easily installed with FWTools. 

http://www.gdal.org/gdal_utilities.html

I've left my partially worked, and "upside-down" example below.

HTH

Cheers, Mike.
To convert with gdalwarp directly from the netcdf file, using the bounding box previously computed in R you will have to reference the subdatasets individually (I'm not how to do this with all the SDS as bands without using VRT virtual driver format): 

         min      max
coords.x1 -8572928 -7012092
coords.x2  4138050  5750943

## first convert the netCDF SDS to the known projection
gdal_translate NETCDF:"20071003.276.0237.n17.nc":mcsst merc-MCSST.tif -a_srs "+proj=merc" -a_ullr -8572928 5750943 -7012092 4138050 

## now convert that to longlat
gdalwarp merc-MCSST.tif longlat-MCSST.tif -t_srs "+proj=longlat"



gdalwarp out.tif Matt Oliver wrote:

Michael, thanks for this example!

One minor question, how would I take the mercator projected ncdf and make it
a longlat projection for the geotiff

(I guess I am wondering at which step I need to do the conversion)

Thanks in advance

Matt

On Thu, Jul 31, 2008 at 7:20 PM, Michael Sumner <[EMAIL PROTECTED]>wrote:


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














_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Reply via email to