Sign up ×
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 files. from the metadata ,the projection is cylindrical and the resolution is 25 km:

                 NC_GLOBAL#institution=SMOS CATDS Processing Chain
                 NC_GLOBAL#references=DDI CATDS - Ref. CAT-DDI-CT-00020-CG
                 NC_GLOBAL#source=CATDS L3 SM Filtrering Processor
             NC_GLOBAL#title=L3-SM P11p : Preliminary Daily Retrieval Map
                               Corner Coordinates:
                    Upper Left  (    0.0,    0.0)
                    Lower Left  (    0.0,  586.0)
                    Upper Right ( 1383.0,    0.0)
                    Lower Right ( 1383.0,  586.0)
                    Center      (  691.5,  293.0)

I want to reproject it to WGS84 with a resolution of 0.25 * 0.25 degree. I used this code:

                    import os
              for doy in range(1,367):
    gdalcmd='gdalwarp -of "ENVI" -ot Float32  -srcnodata -32768 -dstnodata -9999   -te -180 -90 180 90 -overwrite  '+inputDataset+' '+outputFile
     os.system( gdalcmd)

but I got an error:

share|improve this question
The data is probably in the EASE Global grid. The problem is that there are two. One is based upon WGS84 while the other using an authalic sphere based on International 1924 ellipsoid. Can you find any info in the metadata? If it's WGS84, try using 3975 as the input CRS and 4326 as the output. If it's the sphere, try 3410 for the input. –  mkennedy Jan 14 '13 at 19:41
Thanks MKENNEDY ,it is based on authalic sphere based on International 1924 ellipsoid, where to add 3410? –  Jonsson Jan 15 '13 at 14:42
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… has my response. –  Spacedman Jan 19 '13 at 15:56

2 Answers 2

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 - 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

I would suggest to try such conversions with QGIS, giving you a visual impact on what your data looks like. The gdalwarp command basically looks something like

gdalwarp -s_srs EPSG:31466 -t_srs "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" Aerowest.jpg utmtest3.tif

where source and target coordinate reference systems are defined by EPSG codes (the numbers mkennedy provided) or a proj command string.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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