There are a few ways to read the data. Forgive me, but I am not familiar with the java bindings and syntax. One method is to use GDALRasterBand->ReadBlock():
http://gdal.org/classGDALRasterBand.html#a09e1d83971ddff0b43deffd54ef25eef
That reads in the internal block defined by the driver, which is one scanline for the aaigrid driver. See the example in the documentation. To do this manually, similar to your example, use GDALRasterBand->RasterIO(), which appears to be wrapped in java with ReadRaster(). Loop over the raster GDALDataset->GetRasterYSize() times and read in one line:
double[]data = new double[nXSize];
for(int i = 0; i < nYSize; i++)
{
band.ReadRaster(0, i, nXsize, 1, data);
//do something with your data
}
You can change the size of your input window for whatever you would like to get from your data. For example, as above, use a 3x3 window on the ReadRaster() call to calculate slope and aspect. I just found a link to the javadoc here:
http://gdal.org/java/
You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module.
An example:
from math import floor
from osgeo import gdal,ogr
import struct
src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'
src_ds=gdal.Open(src_filename)
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
rb=src_ds.GetRasterBand(1)
ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY() #coord in map units
#Convert from map to pixel coordinates.
px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
px = floor(px) #x pixel
py = floor(py) #y pixel
structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)
print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
Alternatively, since the reason you gave for not using numpy
was to avoid reading the entire array in using ReadAsArray()
, below is an example that uses numpy
and does not read the entire raster in. It uses the built-in gdal.ApplyGeoTransform() function in order to deal with axes rotations.
from math import floor
from osgeo import gdal,ogr
src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'
src_ds=gdal.Open(src_filename)
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
rb=src_ds.GetRasterBand(1)
ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY() #coord in map units
#Convert from map to pixel coordinates.
px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
px = floor(px) #x pixel
py = floor(py) #y pixel
intval=rb.ReadAsArray(px,py,1,1)
print intval[0] #intval is a numpy array, length=1 as we only asked for 1 pixel value
Best Answer
Reading Rasters by block can be done in rasterio and I'd argue it's easier than in GDAL. There is even a tutorial on windowed read/write over at GitHub.
Let's take a look at the
read
functions arguments, which allows you to set a window to read data from:So all you need now is a way to create the window tuples. That's what rasterios
block_windows
function is for:So in your case I guess the code would look something like this: