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