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 for global soil moisture.

From the metadata, the projection is cylindrical and the resolution is 25 km (it is based on authalic sphere based on International 1924 ellipsoid). As I want to compare with other data, I have to make them identical. My other data are in WGS84 with 0.25*0.25 degree resolution and extent -180, 180, -90, 90.

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

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

to:

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 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
    
Spacedman:Thanks ,i have edited the code.where did you post your answer.there is no response from you at all in R-help –  Jonsson Jan 19 '13 at 15:50
    
r.789695.n4.nabble.com/… has my response. –  Spacedman Jan 19 '13 at 15:56
    
Sorry,it is my bad(did not see it). –  Jonsson Jan 19 '13 at 16:01

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

sm=list(x=get.var.ncdf(f,"lon"),y=get.var.ncdf(f,"lat"),z=get.var.ncdf(f,"Soil_Moisture"))
image(sm)

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"){
  require(ncdf)
  require(raster)
  require(sp)
  require(rgdal)

  d = open.ncdf(gridfile)

  sm = list(
    x=get.var.ncdf(d,"lon"),
    y=get.var.ncdf(d,"lat"),
    z=get.var.ncdf(d,name)
    )

  xp = data.frame(x=sm$x,y=0)
  coordinates(xp)=~x+y
  proj4string(xp)=CRS(inCRS)
  xp=spTransform(xp,CRS("+init=epsg:3410"))

  yp = data.frame(x=0,y=sm$y)
  coordinates(yp)=~x+y
  proj4string(yp)=CRS(inCRS)
  yp=spTransform(yp,CRS("+init=epsg:3410"))

  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")
  return(smr)
}

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")
plot(smr)

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,
    nrows=180*4,ncols=360*4,crs="+init=epsg:4326")
  projectRaster(r1,r)
}

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
plot(smrp)

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
    
many thanks for the useful information.Yes the thing is that I want the map to be a regular grid(WGS84). I am grateful to your help and looking forward to further help. –  Jonsson Jan 20 '13 at 10:52
    
Spacedman: I hope you have had time to look at it again –  Jonsson Jan 21 '13 at 13:21
1  
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
1  
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
    
million thanks ,it is deffinitely worth waiting –  Jonsson Jan 22 '13 at 14:59

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