0

Hi I'm working with prism data and I already have clipped (masked) rasters (*.bil files) of a specific area for each day of a year. I'd like to open all these rasters with rasterio, convert them into numpy arrays, get the average value for each array, then an 8 day average of daily arrays for the entire year.

The first step to doing this is to open up all the bil files iteratively, an keep the names of the files as the variable names, but I cannot seem to get this to work.

import glob, os
import fiona
import rasterio
from rasterio.mask import mask
import numpy

path = '/home/demios/Desktop/Prism/Temps'


notebook = dict()

for filename in [each for each in os.listdir('/home/demios/Desktop/Prism/Temps') if each.endswith('.bil')]:
    with open(path+'/'+filename, "r+") as filehandle:
        notebook[filename] = rasterio.open(filehandle)

Near as I can tell this is what my code is supposed start off looking like. It does not seem to do anything though.

  • open with rasterio having matching variable names as files and append variable names to list "rasterload" while opening

  • using last "rasterload convert loaded rasters into numpy arrays with iterative date matching variable names and make a list "arraydates" of them.

  • get average of each "arraydates" into a tuple named "dailyavg" via np.mean(variable)

  • Average dailyavg every 8 days (ie 1-8, 8-16, 16-24, 24-32) into a new tuple 8 day averages.

Anyone able to provide a path forward?

  • Why are you pre-opening the files? You can use rasterio.open() with the path. e.g. rasters = [rasterio.open(p) for p in file_paths]. Then you can stack them with stacked = np.dstack([r.read() for r in rasters]), for instance, assuming their shapes match. – mikewatt Dec 14 '18 at 21:15
  • What do the file names look like? – Marc Pfister Dec 14 '18 at 21:56
  • @MarcPfister they look like "PRISM_tmean_stable_4kmD1_20160101_bil.prj", "PRISM_tmean_stable_4kmD1_20160102_bil.prj", "PRISM_tmean_stable_4kmD1_20160103_bil.prj" etc... Format for values is YYYYMMDD. – Hexadecimalism Dec 17 '18 at 4:21
  • @gberard that is absolutely valid as their shapes match, but I'd still have to find a way to iterate the variable names into some sort of list in order to stack them. – Hexadecimalism Dec 17 '18 at 4:21
  • Ah, I see. If the naming/date convention is consistent then you could just sort your list of filenames and take slices to grab 8 at a time. e.g. for i in range(0, len(images), 8): chunk = images[i:i+8]. Then do what I posted above to stack each set – mikewatt Jan 2 at 23:00
0

Since the file names have the date built into them, it's probably better to build each file name and iterate your way through the year. This way you can build your array in order. Another benefit is that it will error if one of the files is missing.

from datetime import date, timedelta
import calendar
the_year = 2016

# handle leap years
if calendar.isleap(the_year):
    days = 366
else:
    days = 365

# first day of the year you're working on
the_day = date(year=the_year, month=1, day=1)
# a time interval of one day
one_day = timedelta(days=1)

# let's just allocate the whole array since you know the data size in rows and columns
data = np.zeroes((rows, columns, days))

# we want to open one file for every day in 2016:
while the_day.year == the_year:
    # build the path with year/month/day values    
    path = "path/to/PRISM_tmean_stable_4kmD1_{}_bil.prj".format(the_day.strftime("%Y%m%d"))
    # open the file with rasterio
    with rasterio.open(path) as file:
        # the array is zero-indexed, so back up a day
        # .timetuple().tm_yday is the "day of year", so 1-365/6
        index = the_day.timetuple().tm_yday - 1
        #put the data in the right slice of the array
        data[:, :, index] = file.read()
    # advance to the next day
    the_day = the_day + one_day

You can then slice the data array as needed.

  • Sorry for the late comment, but he code only seems to work for the first 31 days (month of january) instead of the whole year, as evidenced by the index maxing out at 30 (0-30 for a total of 31). – Hexadecimalism Feb 21 at 16:28
  • @Hexadecimalism Oh sorry! day.day is the day of month! I have updated it to give the day of the year – Marc Pfister Feb 22 at 1:25
  • Top notch, works like a charm. Point of note "data = np.zeroes(rows, columns, days)" should be "data = np.zeros((rows, columns, days))" – Hexadecimalism Feb 22 at 9:59
  • Ack, I always forget to pass sizes as tuples. Updated with the fix. – Marc Pfister Feb 22 at 17:33

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.