1

I want to use a points shp file to extract some raster files values. It contain two loops: first is raster files loop, second is shp file points loop. The main loop code showed as follows:

allyear={}
for filename in os.listdir(inpath):
    if os.path.splitext(filename)[1]=='.tif':
        ras_file=os.path.join(inpath,filename)

        ras_ds=gdal.Open(ras_file,0)   
        gt=ras_ds.GetGeoTransform()
        rb=ras_ds.GetRasterBand(1)
        arr=rb.ReadAsArray()
        arr=nodata2nan(arr,-32768)   #this is custom function

        print(filename+' is begining processing')

        oneyearvalues={}

        for feat in lyr:
            geom = feat.GetGeometryRef()
            mx,my=geom.GetX(), geom.GetY()  #coord in map units   

            px = int((mx - gt[0]) / gt[1]) #x pixel
            py = int((my - gt[3]) / gt[5]) #y pixel

            pointmeanvalue=bufferpoint(px,py,arr)  #this is custom function
            print(feat['site'],pointmeanvalue,px,py)

            # every feature site is key name and point mean value is key value
            oneyearvalues[str(int(feat['site']))]=pointmeanvalue

         #every raster filename is key name and one extraction is key value
        allyear[os.path.splitext(filename)[0]]=oneyearvalues

total_df=pd.DataFrame.from_dict(allyear)

I checked the printed result,and found that this code just operated the first .tif file. The second for loop just carried out the first raster file,and the others could not be operated. Such as follows:

rasterfile1.tif is begining processing

 ......#processing a short time

rasterfile2.tif is begining processing  
rasterfile3.tif is begining processing
rasterfile4.tif is begining processing
rasterfile5.tif is begining processing

What is the problem?

1

After the first raster is processed all of the features from the layer object have been looped through. For each raster you need to reset the layer object to start at the first feature using the ResetLayer method.

Therefore when you come out of the feature loop, be sure to call lyr.ResetReading():

#...
for feat in lyr:
        geom = feat.GetGeometryRef()
        mx,my=geom.GetX(), geom.GetY()  #coord in map units   

        px = int((mx - gt[0]) / gt[1]) #x pixel
        py = int((my - gt[3]) / gt[5]) #y pixel

        pointmeanvalue=bufferpoint(px,py,arr)  #this is custom function
        print(feat['site'],pointmeanvalue,px,py)

        # every feature site is key name and point mean value is key value
        oneyearvalues[str(int(feat['site']))]=pointmeanvalue

lyr.ResetReading() # Reset the layer back to the first feature

#every raster filename is key name and one extraction is key value
allyear[os.path.splitext(filename)[0]]=oneyearvalues
  • I see! Thx! I know that appearing null value when iterating layer features. – chaobin zhang Jul 13 '17 at 16:04

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.