Writing Numpy array to raster file (tif) returns a trivial black square

Question

I need to turn a 2D numpy array into a raster file that I can later import in ArcGIS (version 10.3.1, licence type: Advanced). The coordinate system that I need to target is the British National Grid coordinates (OSGB1936, EPSG:27700), which are expressed in meters.

This is how I create the tif. This procedure is based on the answer by Eddy The B.

from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import numpy as np
import matplotlib.pylab as plt

#These 3 arrays come from a netCDF file
array=inv_masked
lat=lats #degrees
lon=lons #degrees

xmin,ymin,xmax,ymax = [130859.699132,-964.660232,655859.699132,527035.339768] #meters (obtained with a conversion)
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)

geotransform=(xmin,xres,0,ymax,0, -yres) 

output_raster = gdal.GetDriverByName('GTiff').Create('C:\Users\MyName\Temp\inv_masked.tif',ncols, nrows, 1 ,gdal.GDT_Float32) 
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(27700) #British National Grid OSGB1936
output_raster.SetProjection(srs.ExportToWkt())

The inv_masked array contains boolean values only, and has 352 rows and 350 columns. When plotted with matplotlib, this is what inv_masked looks like:

However, the .tif file that I create with such procedure is just a black square. The spatial reference seems to be correct, but when loaded in ArcGIS I don't see the coastline. Essentially, the .tif map only has 0 values, so all the 1 are not preserved, somehow, and that's why the coastline is not visible. What am I doing wrong here?

What am I doing wrong here?

Answer 1

I cannot see one instruction as this:

output_raster.GetRasterBand(1).WriteArray( array ) 

So, try out next snippet code:

.
.
.

output_raster = gdal.GetDriverByName('GTiff').Create('C:\Users\MyName\Temp\inv_masked.tif',ncols, nrows, 1 ,gdal.GDT_Float32)
#writting output raster
output_raster.GetRasterBand(1).WriteArray( array ) 
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(27700) #British National Grid OSGB1936
output_raster.SetProjection(srs.ExportToWkt())
output_raster = None

Don't forget this line:

output_raster = None

at the bottom of the script.