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

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
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?

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.

• Is there a way to preserve the `(i,j)` indexing of each element of the numpy array when I convert it to `.tif`? Mar 16, 2017 at 13:37
• Index for numpy arrays are inverted as related with GDAL. In GDAL x is related to columns and y is related to rows. Is this the origin of your question? Mar 16, 2017 at 14:25
• Let's see if I can explain myself better. I knew about the `x` and `y` meaning. What I wanted to know is if it is possible, when you save the numpy array as `tif`, to somehow include the `(i,j)` position of each cell in the original numpy array. So that, when you import the raster in ArcGIS or QGIS, you can visualize each cell and retrieve its original position in the numpy array. Mar 16, 2017 at 14:28
• It isn't possible but you can do a small script in QGIS for retrieving its original position in the numpy array with a mouse click on the raster. So, you need to ask another question. Mar 16, 2017 at 14:38
• I have already posted one similar here: gis.stackexchange.com/questions/232350/… Mar 16, 2017 at 14:39