[R-sig-Geo] ncdf file with non-equally spaced grid
Michael Sumner
mdsumner at utas.edu.au
Wed Aug 6 11:53:46 CEST 2008
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 <mdsumner at utas.edu.au>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
>> NC_GLOBAL#creator_email=moliver at udel.edu
>> 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
>> NC_GLOBAL#creator_email=moliver at udel.edu
>> 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 <mdsumner at utas.edu.au
>>
>>> 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: mdsumner at utas.edu.au
>> Phone: 03 6226 1752 (W)
>> 0408599921 (M)
>> Fax: 03 6226 2745
>>
>>
>>
>>
>>
>>
>>
>>
>>
>
>
>
More information about the R-sig-Geo
mailing list