# Generating a raster file using 1D numpy arrays

Explanation: I currently have 3 1D numpy arrays, each having a length of over 100000 points. The first contains the X coordinates of an object, the second contains the Y coordinates, and the third one contains the speed of the object at the corresponding coordinates.The 3 arrays are in sequence from the start to end. I am trying to generate a raster file that will match another raster file that I produced using certain images.For simplicity's sake, the output raster size is 400 meters x 400 meters, and each cell has dimensions of 20 X 2. How exactly do I go about generating a raster using my 3 arrays ?

Ideally I would like to do something as simple as:

``````#xm = x-cooridnate array
#ym = Y-coordinate array
#speed = speed array

xmin,ymin,xmax,ymax = [xm.min(),ym.min(),xm.max(),ym.max()]
#nrows,ncols = np.shape(mining_durations)
xres = (xmax-xmin)/float(2)
yres = (ymax-ymin)/float(20)
geotransform=(xmin,xres,0,ymax,0, -yres)
output_raster = gdal.GetDriverByName('GTiff').Create(outputFileName,400/20, 400/2, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)
output_raster.GetRasterBand(1).WriteArray(speed)
output_raster.FlushCache()
``````

But this of course doesn't work.I know that many points will coincide in one cell, but I assume that it will automatically take the average of all the coinciding points.

Do I need to first create some sort of grid for the speed array?

Edit: I have added sample data, as well as the code I am using.

Sample data: The coordinates are just float values. For example, the first column here is X, the second is Y and the last is some other variable (speed):

``````52039,57 -54169,66 7,68667678522605
52027,68 -54180,01 6,8999999999869
52039,55 -54169,70 2,12132034354421
52039,59 -54169,71 1,34164078652916
52039,61 -54169,64 0,670820393288984
52039,63 -54169,51 0,600000000013097
52039,69 -54169,49 0,750000000043656
52039,70 -54169,46 0,150000000030559
52039,67 -54169,58 0,75
``````

Each sample sample/row of data is taken every 4 seconds, although this sometimes changes. The X and Y coordinates don't necessarily follow a fixed pattern. The code I am using is shown below.Note that this doesn't produce an output of 400x400 which I originally had in mind, the size is based on the coordinates.

``````def GenerateRaster(xdata, ydata, zdata):

binWidth = 2.0
binLength = 20.0
xMin = min(xdata)
xMax = max(xdata)
yMin = min(ydata)
yMax = max(ydata)

statistic,x_edge,y_edge,binnumber = stats.binned_statistic_2d(xdata, ydata, zdata, 'mean',bins = [np.arange(xMin, xMax, binWidth), np.arange(yMin, yMax, binLength)] )
nrows,ncols = np.shape(statistic)

xres = (xMax-xMin)/y_edge.shape
yres = (yMax-yMin)/x_edge.shape

geotransform=(xMin,xres,0,yMax,0, -yres)
output_raster = gdal.GetDriverByName('GTiff').Create('Test_Raster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
output_raster.GetRasterBand(1).WriteArray(statistic)   # Writes array to the raster
output_raster.FlushCache()
``````

Although I do produce a raster file, I am not sure if this is correct or accurate. I noticed that the cell sizes are not exactly 20 x 2 which is strange: the dimensions are slightly longer.

• Can you explain the data in the three arrays in a little more detail (and perhaps include a few rows)? What coordinate reference system are the coordinates in? Are the points equally spaced? Do the points fully cover the 400 m x 400 m area? Are there any points outside the 400 m x 400 m area? What are the coordinates for the corners of the 400 m x 400 m area? Jun 12, 2018 at 4:42
• Hi Dave. I was not given a coordinate reference system so when I open the tiffs in gdal I just use the default one that comes up. I have updated the post to include a sample of data and the code I am currently using to produce the raster. This isn't in the 400m x 400m area I originally wanted, instead its based on the minimum and maximum values, which is what it should be, right ? Even if it doesn't match my other rasters perfectly, the X and Y positions at a particular point should give an approximate comparison, at least I hope it does. Jun 12, 2018 at 5:59
• It seems you answered your initial question regarding the need for a grid with your use of stats.binned_statistic_2d function for that purpose. I think your cell size anomoly is due to 2 things: (1) `np.arange(xMin, xMax, binWidth)` will produce `[xMin, xMin+binWidth, ..., xMax - binWidth]` (i.e., your bins do not fully cover the data range). Try `np.arange(xMin, xMax+binWidth, binWidth)`. And (2) for the resolution calcs: I think you want the number of bins rather than the number of edges in the denominator. Also, shouldn't `xres` have `**x**_edge.shape` in the denominator? Jun 12, 2018 at 17:31
• Your first suggestion has partially worked, thank you! With regards to your second suggestion: by "number of bins" do you mean number of columns or rows ? So `xres = (xMax-xMin)/ncols`?. What do you mean by "`**x**`" ? I think `xres` has to have a `y_edge.shape` or `ncols` in the denominator, because if has i had `xres = (xMax-xMin)/x_edge.shape` or `yres = ../y_edge.shape` then it would be stretched. After playing around with different resolutions, I think I am going to add your first suggestion, as the cell dimensions of the raster were very close to 20 x 2 (although not exactly). Jun 13, 2018 at 7:33
• Sorry, `**x**` was a formatting goof. By bins, I did mean ncols & nrows - however, that won't work either. Consider one axis: first bin opens on the min value and subsequent bins open and close on an interval per the bin size. It's unlikely that the max value is exactly a multiple of the bin size + the min value. As such, the last bin closes past the max value. Therefore, a resolution based on min & max would be close, but likely wrong. Should use xmin & the close of the last bin. But resolution is bin size! TLDR: use `binWidth` & `binLength` as the resolutions for creating the tiff Jun 13, 2018 at 21:56

`np.arange(xMin, xMax, binWidth)` will produce `[xMin, xMin+binWidth, ..., xMax - binWidth]` (i.e., your bins do not fully cover the data range).
Try `np.arange(xMin, xMax+binWidth, binWidth)`