rsgislib.zonalstats.pixelVals2TXT()

[Sam Alaee]

Hello,

I am looking for a way to get all the pixel values of each of the polygons for all of the bands. The closest thing that I could find is using (rsgislib.zonalstats.pixelVals2TXT()) but this function exports csv files for each polygon and it is computationally expensive when there are large number of segments, and it takes a lot of time. I need to be able to select a certain polygon and a specific band in my vector file and get the pixel values in form of an array. 




I was wondering if any of you could help me with this. 




[Nathan Thomas]

Hi Sam,

rsgislib.imageutile.extractZoneImageValues2HDF will be faster than pixelVals2TXT. It will require that your polygons are rasterized to the same resolution as your input image as a binary image (1 or 0). It will extract values to a hdf5 file. You can then use the python module h5py to read the data from the hdf5 file directly as a numpy array.

steps:

1. rasterize polygons with gal_rasterize (or QGIS)

2. run rsgislib.imageutils.extractZoneImageValues2HDF(‘Landsat_sref.kea’, ‘RasterizedPolygon1.kea’, 'ExtractionPolygon1.h5’, 1)

3. open values in python:

import h5py

f = h5py.File( 'ExtractionPolygon1.h5’, ‘r’)

dset = f[‘DATA/DATA’]

band1 = dset[…,0]

band2 = dset[…,1]

Nathan

Nathan M Thomas

CalTech Postdoc at JPL

--
You received this message because you are subscribed to the Google Groups "RSGISLib Support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-suppo...@googlegroups.com.
To post to this group, send email to rsgislib...@googlegroups.com.
Visit this group at https://groups.google.com/group/rsgislib-support.
For more options, visit https://groups.google.com/d/optout.

[Sam Alaee]

Hi Nathan,

Thanks. I tried your approach, it's much better than pixelVals2TXT. I was wondering if there was a way that didn't need to create any files. 

[Nathan Thomas]

Hi Sam,

Im not aware of a way within RSGISlib to do it without creating additional files. Pete or Dan may know of a way, but there isn’t one to my knowledge. You could read in the image and the rasterized polygon using gdal, but even then you would have to make sure that the two datasets are on exactly the same grid, which would create additional files.

You would only have to do this once so you could try:

1. rasterize polygons to same grid as imagery (specify the image extent with -te)

2.

        from osgeo import gdal

    driver = gdal.GetDriverByName('KEA')

    image = ‘landsat.kea'

    poly = ‘polygon.kea'

    

    Idataset = gdal.Open(image)

    Pdataset = gdal.Open(poly)

    # Landsat band 1. Needs to changed for other bands

    Iband = Idataset.GetRasterBand(1)

    # polygon only has band 1

    Pband = Pdataset.GetRasterBand(1)

    # Get Rows and Cols numbers for each image

    ICols = Idataset.RasterXSize

    IRows = Idataset.RasterYSize

    PCols = Pdataset.RasterXSize

    PRows = Pdataset.RasterYSize

    # Read images in as an array

    print('Reading Radar data...')

    Idata = Iband.ReadAsArray(0,0,ICols,IRows).astype(numpy.float)

    print('Reading index data...')

    Pdata = Pband.ReadAsArray(0,0,PCols,PRows).astype(numpy.float)

    

    # subset landsat image where polygon values are valid (i.e 1)

    values = Idata[Pdata==1]

Nathan

    

Nathan M Thomas

CalTech Postdoc at JPL

Radar Science and Engineering

4800 Oak Grove Dr.

300-356A

Pasadena, CA 91109

Phone: +1 818 354 4681

[Sam Alaee]

Hi Nathan,

Awesome. Thanks for the help.

Best,

Sam

To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-support+unsub...@googlegroups.com.

[Daniel Clewley]

If you don’t want to keep the additional files you could always use the tempfile module (https://docs.python.org/3/library/tempfile.html) to create a temporary file / directory and them remove it at the end of the script (e.g., `os.remove(tempfile)`).

Thanks,

Dan