1

Is there a very simple way to plot up variables from a netCDF file that are stored in a regular lat-lon grid (-90to90˚N, 0-360˚E) onto the now more popular Mollweide projection?

For example say I have something like the following:

lon <- seq(-90,90,2.5) # really these data would come from the netCDF file...
lat <- seq(0,356.25,3.75)
myvar <- runif(length(lon)*length(lat))
mymat <- array(dim=c(length(lon), length(lat)), myvar)

Which I can plot up easily with something like image:

image(lon, lat, myvar)

But what I'd like to do is convert the myvar array and reproject it so that I can plot it in the Mollweide projection:

... convert to Mollweide projection ...
image(myvar_converted)

I'm looking for a very generic approach and ideally something that would allow me to add on points, lines etc on to the plots using base functions if possible.

  • You should start giving a look to the following R packages: raster, rgdal and maps. Then you'll realize that, thanks to Google, there are many tutorials available, like this: blog.fellstat.com/?p=130 – Matteo De Felice Aug 14 '15 at 9:37
  • Take a look at the gdalUtils package, specifically gdalwarp. You'll need to write your source data out first though. – MikeRSpencer Aug 16 '15 at 21:21
1

I'll just post a few pointers here for now: when I get a chance, I'll try to make this a more comprehensive and easily-consumable answer. But for now:

First off, consider changing your title better match (IIUC--ICBW) the terms folks in this domain use for this sort of usecase, and therefore to make it more findable. What you want to do is more usually called projecting spatial data from unprojected/lon-lat (aka "global") to Mollweide (or other desired output projection). (This may also be called "reprojecting" or "regridding," but "projecting" is more correct for this usecase.)

I have code (which needs refactoring, but works and is IMHO rather better documented than a lot of science code on the web) that does something similar to what you want (along with lots of other things that you don't want): projecting netCDF from global lon-lat to a different projection (in my case, LCC), e.g.:

Of these, GEIA_regrid is the oldest, and probably least well-written, but it does use (IIRC) R package=image, which you seem to want. My recommendations:

  1. Read up on the different R plot packages. My experience was, I started with image but rapidly found I needed lattice functionality, and eventually realized that I should have started by learning ggplot2. There now may be an even newer/better R package for spatial plotting. (If your project is truly a one-off this won't matter.)
  2. Read the READMEs for either GFED-3.1_global_to_AQMEII-NA (which IIRC is a bit newer/cleaner, but a bit more complex due to the way GFED provides its data) or CLM_CN_global_to_AQMEII-NA to (hopefully :-) get some idea how these codes/projects work. To a first approximation, these projects
    1. use bash to drive NCL and R code. (If you're gonna be working with netCDF much, consider learning a bit about NCL: it's "the bomb" for netCDF hacking.) If I was starting over, I might try an all-python solution, since the python netCDF libraries are pretty good.
    2. have an uber_driver.sh which sets up the environment and filespace (either from local filetree or from git repo) and then calls a "project driver" (e.g., CLMCN_driver.sh, GFED_driver.sh), which calls the NCL and R work scripts.
  3. The work scripts that most do what you want are the *regrid* ones: regrid_global_to_AQMEII.r, vis_regrid_vis.r. You'll probably need to make a buncha changes to make these do what you want, but you will definitely need to change one thing: output-projection specification in variable name=out.proj4. This is a PROJ.4 string, which is "the standard way" to really, really specify geographic projections, or lack thereof (in the lon-lat case).
    • I was required to work with a rather "fiddly" output format, so I set out.proj4 programmatically using a library that deeply understands my output format.
    • You will probably want to explicitly set your out.proj4 for Mollweide as indicated here, but (as indicated there) you will need to know the appropriate false easting, false northing, and central meridian for your output in order to set that.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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