1

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: enter image description here

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?

enter image description here

What am I doing wrong here?

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.

| improve this answer | |
  • Is there a way to preserve the (i,j) indexing of each element of the numpy array when I convert it to .tif? – FaCoffee Mar 16 '17 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? – xunilk Mar 16 '17 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. – FaCoffee Mar 16 '17 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. – xunilk Mar 16 '17 at 14:38
  • I have already posted one similar here: gis.stackexchange.com/questions/232350/… – FaCoffee Mar 16 '17 at 14:39

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.