Friday 14 April 2017

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



Answer



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.



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...