1

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)] )
   # statistic = np.ma.masked_invalid(statistic)
    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? – Dave G Jun 12 '18 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. – user3451660 Jun 12 '18 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? – Dave G Jun 12 '18 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). – user3451660 Jun 13 '18 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 – Dave G Jun 13 '18 at 21:56
0

The useful bits from my comments are below.

Regarding the cell sizes:

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.

Regarding portion of data being missed:

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.