Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have nc (NetCDF) files .

From the metadata, the projection is cylindrical and the resolution is 25 km (it is based on authalic sphere based on International 1924 ellipsoid).

I want to re-project the data I have from:

EASE grid, cylindrical, 25 km, 1383 pixel 586 lines


WGS84, 0.25*0.25 degree, 1440 pixel 720 lines

share|improve this question

migrated from stackoverflow.com Jan 19 '13 at 14:17

This question came from our site for professional and enthusiast programmers.

marked as duplicate by whuber Apr 1 '14 at 16:16

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

Downvote for: a) not reproducible (is the .nc file the same as the .dbl file in the download?) b) cross-posting to R-help c) ignoring my response in R-help. d) not telling us package version numbers –  Spacedman Jan 19 '13 at 15:34
r.789695.n4.nabble.com/… has my response. –  Spacedman Jan 19 '13 at 15:56

1 Answer 1

up vote 5 down vote accepted

I have a vague outline of how to do this, but there's plenty that you may have to understand.

NetCDF files are complex general data containers so its not always clear how to get spatial data out of them. In this case, you can get the Soil_Moisture variable and that is just a 2d matrix with no coordinate reference. If you do image(A) you should see your soil moisture data, but the X and Y axes won't be correct.

Your particular NetCDF files have lon and lat members, which are the coordinates of the cells. These are NOT in EASE coordinates, they are in lat-long coordinates. If you do


Then you will see your data with correct lat-long coordinates. This overlays quite nicely with lat-long data from things in WGS84 coordinates.

BUT... your grid is not composed of square pixels. Consecutive y-coordinates do not have a constant difference. This is because they define pixels of different shape as you go to the poles. This means that you can't just create a raster:

> smr=raster(sm)
Error in .local(x, ...) : data are not on a regular grid

which is the first step to transformation.

What you need to do is to convert those x and y coordinates to values in the EASE coordinate system - which seems to be EPSG code 3410 http://www.spatialreference.org/ref/epsg/3410/ - and then they should be a regular grid with a constant step between them - the units just won't be degrees or metres but projected units.

Here's that code. I'm not sure what the lat-long coordinates given in the NetCDF file is. It could be WGS84 or it could be something else. The other likely candidate looks like EPSG code 4053, but they both give the same answer. At this scale, the exact earth shape probably doesn't matter that much. I've made that CRS a function so you can try other ones easy enough.

Note how I get all the grid X values and create a dataframe with zeroes and convert that to get all the X values in the EPSG 3410 system, then do a similar thing with the Y values. That should get a regular spaced numeric system that the raster function can handle:

convertGrid <- function(gridfile, name, inCRS="+init=epsg:4053"){

  d = open.ncdf(gridfile)

  sm = list(

  xp = data.frame(x=sm$x,y=0)

  yp = data.frame(x=0,y=sm$y)

  sm$xp = coordinates(xp)[,1]
  sm$yp = coordinates(yp)[,2]

  smr = raster(list(x=sm$xp,y=sm$yp,z=sm$z),crs="+init=epsg:3410")

So then I can read in the soil moisture as a raster:

smr = convertGrid("SM_RE01_MIR_CLF31D_20100812T000000_20100812T235959_246_001_7.DBL","Soil_Moisture")

Now to get onto the coordinate system you really want. First define an empty raster with the desired number of cells and location, then use projectRaster to interpolate the cells to the new basis.

transformTo <- function(r1){
### 0.25*0.25 degree resolution and extent -180, 180, -90, 90
  r=raster(xmn=-180, xmx=180, ymn=-90, ymx=90,

Now use that function to transform, first setting missing values:

## anything below zero is NA (-1 is missing data, soil moisture is +ve)
smr[smr < -0.1] <- NA
smrp = transformTo(smr) # takes a short while

Now smrp should be in the exact grid that you want.

> smrp
class       : RasterLayer 
dimensions  : 720, 1440, 1036800  (nrow, ncol, ncell)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 1  (min, max)

Worth waiting for?

share|improve this answer
Spacedman: I hope you have had time to look at it again –  Jonsson Jan 21 '13 at 13:21
Quick note though: R seems to read Soil_Moisture with values -1 to 1, but the ncview program shows values above 0 with -1 marking missing data. Ah, I think the interpolation has interpolated between the 0 and the -1. You need to set anything below 0 in the raster to NA... –  Spacedman Jan 21 '13 at 22:36
Also I think the reason epsg:4326 for lat-long doesn't change anything is because Proj4 isn't doing a datum shift here. I wouldn't worry too much... –  Spacedman Jan 21 '13 at 22:50

Not the answer you're looking for? Browse other questions tagged or ask your own question.