Hi Kate,
There are a couple of ways you can skin this cat. An easy and computationally efficient way to do it without the Mapping Toolbox is just treat longitudes like x values and latitudes like y values. Then use my borders function to draw national boundaries:
imagesc(lon,lat,cummDose)
borders('countries','nomappingtoolbox')
axis xy
will plot everything you want to plot in plain unprojected coordinates. Since you seem to have some missing data, you may wish to make it transparent with imagesc by
h = imagesc(lon,lat,cummDose);
set(h,'alphadata',isfinite(cummDose))
If you prefer to do this with the Mapping Toolbox, you'll lose some computational efficiency and it'll be harder to share your scripts with anyone who does not have the Mapping Toolbox, but you'll gain the nice map projections available via worldmap.
Your dataset is not huge, so computational efficiency won't be a big issue, but if you ever experience figure sluggishness, try doing it with imagesc instead of worldmap. The difference is imagesc plots an easy, predictable grid, whereas projecting to map coordinates or using pcolor lets the grid squeeze and stretch. That said, to do this with the Mapping Toolbox, start by initializing a map with worldmap, then use pcolorm, then add the landareas as follows:
worldmap('world')
pcolorm(lat,lon,cummDose)
geoshow('landareas.shp')
