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())
.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?