4

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?

Community
  • 1
FaCoffee
  • 823
  • 8
  • 21

1 Answers1

3

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.

xunilk
  • 29,891
  • 4
  • 41
  • 80
  • 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: http://gis.stackexchange.com/questions/232350/arcgis-raster-table-only-shows-two-groups-of-cells – FaCoffee Mar 16 '17 at 14:39
  • But the question is for ArcGis. However, it is possible for QGIS; as I said above. Could you ask in relation to QGIS? – xunilk Mar 16 '17 at 14:51
  • All right, I will modify that so to accommodate both solutions. – FaCoffee Mar 16 '17 at 14:52
  • Remember that the question can be closed if it has more than one subject. – xunilk Mar 16 '17 at 15:01