python - Turn a shapefile into a mask and calculate the mean?

I have been using Python to process huge arrays stored as NetCDF files. I would like to calculate the mean of an area defined by a shapefile. I have just installed GDAL but if there are other tools I should use please let me know.

So my main question is how can I turn my shapefile into a mask?

Then use my mask to calculate the NumPy mean of my array within the area of the mask?

shp="country.shp"
array=mynumpyarray
mask=useGDALtochangeshapefiletomask
meanofarea=N.mean(array,limits=mask)

Edit

I have been able to convert my shapefile to a raster using gdal_RasterizeLayer link supplied and then a masked array using

maskarray=mask_ds.GetRasterBand(1).ReadAsArray()

Now it is just a question of matching my NetCDF extent with my raster/maskarray. Here are my extents as requested.

my shapefile/raster extent: x_min, x_max, y_min, y_max 140.962408758 149.974994992 -39.1366533667 -33.9813898583

my netcdf file extent: min longitude: 139.8 max longitude: 150.0 min latitude -39.2 max latitude: -33.6 LAT size 106 LON size 193