Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have approximately 700 raster files in a single folder and I need to find files that have the same spatial extent and copy to a different folder. The raster files measure daily snow cover and therefore have the exact same spatial extent.

share|improve this question
1  
Do you have any code assembled? Any particular package you want to use, such as ArcGIS/arcpy, GDAL, etc.? –  Chad Cooper Dec 9 '11 at 13:41
    
I don't have code assembled. I'm not sure where to start as I'm fairly new to using python. –  dchaboya Dec 9 '11 at 15:41
add comment

1 Answer

One possible solution using GDAL.

The first trick is get a list of rasters into Python. For argument's sake we'll assume they are all TIFF files, meaning we can use glob to parse the directory:

import glob
raster_list = glob.glob("path_to_folder/*.tif")

I'm going to assume that you don't know what the size is, so you want different folders for each set of rasters with the same size. We'll use the GDAL's RasterXSize, RasterYSize dataset attributes and the GetGeoTransform() method to compare Rasters (see Week 4 slides here for examples of how these work), storing the results in a dictionary.

import shutil
import os
from osgeo import gdal
gdal.AllRegister() #Registers data formats

i = 0 #We'll start with folder zero
extents = {}

for r in raster_list:
    ds = gdal.Open(r, GA_ReadOnly) #GA_ReadOnly is the GDAL Access level
    gt = ds.GetGeoTransform()
    top_left_x = gt[0]
    top_left_y = gt[3]
    bottom_right_x = top_left_x + ds.RasterXSize * gt[1] #x cell size
    bottom_right_y = top_left_y - ds.RasterYSize * gt[4] #y cell size
    #we use a tuple of the extent to store the destination folder
    if (top_left_x, top_left_y, bottom_right_x, bottom_right_y) not in extents:
        f_name = "extent" + str(i) #new folder name
        i += 1
        os.mkdir(f_name)
        extents[(top_left_x, top_left_y, bottom_right_x, bottom_right_y)] = f_name
    #regardless, we copy the file
    shutil.copyfile(r, extents[(top_left_x, top_left_y, bottom_right_x, bottom_right_y)])

Note you could perform an almost identical operation with arcpy.ListDatasets("*", "Raster") and arcpy.Describe(raster).extent if you have ArcGIS (see here and here for docs).

share|improve this answer
    
Hmmm, slight misswording on what I was asking (6 a.m. posts are not recommended). What I have are rasters that are immediatly adjacent to one another (think of it as a square that was divided into four smaller squares.) now I need to find all rasters with the same spatial coordinates. Extent will not work as all 700 rasters have the same size. Is it possible to generate four polygons that cover each of the four different coordinates and search for rasters using those polygons? –  dchaboya Dec 9 '11 at 15:50
    
I can see why just number of cells wouldn't work if the rasters are the same size, but shouldn't you get four different sets of extents based on which quarter the raster is in (in projection units)? I could be reading this wrong though. –  om_henners Dec 12 '11 at 9:56
    
Good idea, I didn't of using the projection units for the extent parameters. Thanks for the suggestion. –  dchaboya Dec 12 '11 at 13:07
    
No worries. If it helps, the example above the x origin (top_left_x) and the x cell size are in projection units. The RasterXSize is the actual cell count. Likewise for y. So in theory it does use the projected extent to calculate overlap. Good luck! –  om_henners Dec 12 '11 at 23:53
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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