Vjylku

Exploding Numbers

Given the mapping a = 1, b = 2, ... z = 26, and an encoded message, count the number of ways it can be decoded

How Create a list of the first 10,000 digits of Pi and sum it?

Why does this quiz question say that protons and electrons do not combine to form neutrons?

Does this property characterize the odd primes?

Is it possible to detect 100% of SQLi with a simple regex?

What does it mean when an external ID field follows a DML Statement?

What is an explicit bijection in combinatorics?

Are encryption algorithms with fixed-point free permutations inherently flawed?

How can guns be countered by melee combat without raw-ability or exceptional explanations?

Why don't you get burned by the wood benches in a sauna?

Are Kickstages exclusive to Rocketlabs Electron Design?

How do I avoid the "chosen hero" feeling?

How does holding onto an active but un-used credit card affect your ability to get a loan?

Can a Hydra make multiple opportunity attacks at once?

Did the characters in Moving Pictures not know about cameras like Twoflower's?

How many copper coins fit inside a cubic foot?

Can I do anything else with aspersions other than cast them?

Isn't a semicolon (';') needed after a function declaration in C++?

Can I legally make a website about boycotting a certain company?

What does an unprocessed RAW file look like?

Build ASCII Podiums

Will linear voltage regulator step up current?

Badly designed reimbursement form. What does that say about the company?

Loading all rasters in a directory and naming the variables iteratively

Writing numpy array to raster fileUsing GDAL/Python to stack georeferenced images of different sizesFast raster combine using rasterioCalculating Climate Suitable Days using netCDF4 and Python?Build a table with columns of raster values from multiple raster datasets, using Python, GDAL, or PyQGIS?Rasterio + QGIS - rasterio.drivers() -> TypeError: 'module' object is not callableEfficient way to analyze multiple (large) rasters against a different raster using numpyArcGIS 10.2 Raster Calculator can't find numpy after new local Python was installed separatelyComputing mean of all rasters in a directory using pythonConverting ESRI Grid to GeoTIFF and writing to NetCDF

0

1

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?

python raster numpy rasterio

share|improve this question

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  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
  

  
  
  
  

add a comment | 

0

1

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?

python raster numpy rasterio

share|improve this question

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  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
  

  
  
  
  

add a comment | 

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?

python raster numpy rasterio

share|improve this question

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

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)

share|improve this question

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

share|improve this question

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

asked Dec 14 '18 at 19:12

HexadecimalismHexadecimalism

54

1 Answer
1

active

oldest

votes

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

the_year = 2016



# 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, 365)



# we want to open one file for every day in 2016:

while the_day.year == the_year:

    # build the path with yea/mont/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

        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.

share|improve this answer

edited 7 mins ago

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  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
  9 hours ago
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @Hexadecimalism Oh sorry! day.day is the day of month! I have updated it to give the day of the year
  
  – Marc Pfister
  7 mins ago
  

  
  
  
  

add a comment | 

Your Answer

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "79"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});


}
});

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fgis.stackexchange.com%2fquestions%2f306098%2floading-all-rasters-in-a-directory-and-naming-the-variables-iteratively%23new-answer', 'question_page');
}
);

Post as a guest

Name

Email

Required, but never shown

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

the_year = 2016



# 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, 365)



# we want to open one file for every day in 2016:

while the_day.year == the_year:

    # build the path with yea/mont/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

        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.

share|improve this answer

edited 7 mins ago

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  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
  9 hours ago
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @Hexadecimalism Oh sorry! day.day is the day of month! I have updated it to give the day of the year
  
  – Marc Pfister
  7 mins ago
  

  
  
  
  

add a comment | 

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

the_year = 2016



# 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, 365)



# we want to open one file for every day in 2016:

while the_day.year == the_year:

    # build the path with yea/mont/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

        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.

share|improve this answer

edited 7 mins ago

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  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
  9 hours ago
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @Hexadecimalism Oh sorry! day.day is the day of month! I have updated it to give the day of the year
  
  – Marc Pfister
  7 mins ago
  

  
  
  
  

add a comment | 

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

the_year = 2016



# 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, 365)



# we want to open one file for every day in 2016:

while the_day.year == the_year:

    # build the path with yea/mont/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

        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.

share|improve this answer

edited 7 mins ago

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799

from datetime import date, timedelta

the_year = 2016



# 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, 365)



# we want to open one file for every day in 2016:

while the_day.year == the_year:

    # build the path with yea/mont/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

        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

share|improve this answer

edited 7 mins ago

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799

answered Dec 17 '18 at 15:42

Marc PfisterMarc Pfister

3,22799