python - How to re-project the EASE (Equal Area Scalable Earth) grid with a ~25 km cylindrical projection to WGS84 0.25 degree?

I have nc files. from the metadata ,the projection is cylindrical and the resolution is 25 km:

           indrical
                 NC_GLOBAL#ease_resolution=25
                 NC_GLOBAL#history=none
                 NC_GLOBAL#institution=SMOS CATDS Processing Chain
                 NC_GLOBAL#netcdf_version_id=3.6.2
                 NC_GLOBAL#product_version=1.0
                 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):
inputFile='D:\\358\\2010SMOSl3_sm'+str(doy)+'.nc'
inputDataset='NETCDF:\"'+inputFile+'\":Soil_Moisture'
    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:

Answer

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?