'''The following script samples NETCDF files across time and space to sample temperature and preciptation data''' #---------------------------# #-------- LIBRARIES---------# #---------------------------# from osgeo import gdal from osgeo.gdalnumeric import * from osgeo.gdalconst import * import numpy import math import random import glob import csv from itertools import izip_longest from osgeo.osr import SpatialReference, CoordinateTransformation #---------------------------# #---------FUNCTIONS---------# #---------------------------# def make_lat_long_list(fileName,return_column_name): # open the file in universal line ending mode with open(fileName, 'rU') as infile: reader = csv.DictReader(infile) data = {} for row in reader: for header, value in row.items(): try: data[header].append(value) except KeyError: data[header] = [value] return data[return_column_name] def concatenate_lists(column1,column2): 'a simple concatenate function' final_list = [] for i in range(len(column1)): final_list.append(column1[i]+str(' ')+column2[i]) return final_list def sample_tif(tif,latitude,longitude): 'samples a tif based on latitude and longitude and writes the file' source_file = tif src_ds=gdal.Open(source_file) gt = src_ds.GetGeoTransform() rb=src_ds.GetRasterBand(1) temperatures = [] for i in range(len(latitude)): plong = int((180.0 - abs(float(longitude[i])))/ gt[1]) plat = int((90.0 - float(latitude[i]))/ gt[1]) #y interval = rb.ReadAsArray(plong,plat,1,1) temperatures.append(interval[0][0]) with open('temperatures.csv','w') as f: writer = csv.writer(f); writer.writerows(zip(temperatures)); def generate_monthly_mins_maxs(directory,latitude,longitude,species_names,masses): 'A function to generate monthly minimum and maximum temperatures from 30 year intervals' maps = glob.glob('Documents/manuscripts/grinnell_resurvey_project/map/climate_data/grinnell_worldclim/all/*tif') list_of_months = [] for each_month in range(len(maps)): source_file = maps[each_month] src_ds=gdal.Open(source_file) gt = src_ds.GetGeoTransform() rb=src_ds.GetRasterBand(1) temperatures = [] for i in range(len(latitude)): plong = int((180.0 - abs(float(longitude[i])))/ gt[1]) plat = int((90.0 - float(latitude[i]))/ gt[1]) interval = rb.ReadAsArray(plong,plat,1,1) temperatures.append(interval[0][0]) list_of_months.append(temperatures) zip_T = list(izip_longest(*[species_names,masses,latitude,longitude,list_of_months[0], list_of_months[1], list_of_months[2], list_of_months[3], list_of_months[4], list_of_months[5], list_of_months[6], list_of_months[7], list_of_months[8], list_of_months[9], list_of_months[10], list_of_months[11], list_of_months[12], list_of_months[13], list_of_months[14], list_of_months[15], list_of_months[16], list_of_months[17], list_of_months[18], list_of_months[19], list_of_months[20], list_of_months[21], list_of_months[22], list_of_months[23]])) with open('temp_min_max_db.csv','w') as f: writer = csv.writer(f); writer.writerow(['species','mass','latitude','longitude','mo1_max','mo2_max','mo3_max','mo4_max','mo5_max','mo6_max','mo7_max','mo8_max','mo9_max','mo10_max','mo11_max','mo12_max','mo1_min','mo2_min','mo3_min','mo4_min','mo5_min','mo6_min','mo7_min','mo8_min','mo9_min','mo10_min','mo11_min','mo12_min']) writer.writerows(zip_T); def sample_netCDF(directory,latitude,longitude,year): 'a function that sample NetCDF files by year' map_ = glob.glob(directory+str('*.nc')) source_file = map_[0] src_ds=gdal.Open(source_file) gt = src_ds.GetGeoTransform() #gt is the spatial resolution annual_temperatures = [] for each_point in range(len(latitude)): monthly_temperatures = [] firstmonth = (float(year[each_point]) - 1900)*12 for i in range(12): index = int(1 + firstmonth + i) rb=src_ds.GetRasterBand(index) #rb is the specific raster band plong = int((360.0 - abs(float(longitude[each_point])))/ gt[1]) plat = int((90.0 - float(latitude[each_point]))/ gt[1]) interval = rb.ReadAsArray(plong,plat,1,1) monthly_temperatures.append(interval[0][0]) if each_point == 0 or each_point == 29000: print(monthly_temperatures) annual_temperatures.append(mean(monthly_temperatures)) zip_T = list(izip_longest(*[latitude,longitude,year,annual_temperatures])) with open('temperatures.csv','w') as f: writer = csv.writer(f); writer.writerow(['latitude','longitude','year','temperature']) writer.writerows(zip_T); def sample_netCDF_CA(longitude,latitude): 'a function that sample NetCDF files by year for the state of California'\ 'key pint: adjust the years (modern or historic) and name of files individually' #(-374495.83635354, 270.0, 0.0, 592636.6658113, 0.0, -270.0) map_ = glob.glob('*.nc') source_file = map_[0] src_ds=gdal.Open(source_file) gt = src_ds.GetGeoTransform() #gt is the spatial resolution #years = range(0,30,1) years = range(0,115,1) #years = range(85,115,1) months = range(1,13,1) coords_list = [] for each_coord in range(len(latitude)): yearly_list = [] for month in range(len(months)): monthly_value = 0 for every_year in range(len(years)): index = (years[every_year]*12) + months[month] lat = float(latitude[each_coord]) lon = float(longitude[each_coord]) rb=src_ds.GetRasterBand(index) #rb is the specific raster band plat = int((gt[3] + abs(lat))/gt[1]) plong = int((lon + abs(gt[0]))/gt[1]) if plat > 4477 or plong > 3486: pass else: interval = rb.ReadAsArray(plong,plat,1,1) value = interval[0][0] monthly_value += value/100.0 yearly_list.append(monthly_value/len(years))#this is where you calculate the average coords_list.append(yearly_list) zip_T = list(izip_longest(*[latitude,longitude,coords_list])) #with open('CA_Temps_1896_1925.csv','w') as f: #with open('CA_Temps_1996_2015.csv','w') as f: with open('CA_Temps_1896_2015.csv','w') as f: writer = csv.writer(f) writer.writerow(['latitude','longitude','temperature']) writer.writerows(zip_T) def sample_netCDF_CA_by_year(latitude,longitude): 'a function that sample NetCDF files by year for the state of California'\ #(-374495.83635354, 270.0, 0.0, 592636.6658113, 0.0, -270.0) map_ = glob.glob('*.nc') source_file = map_[0] src_ds=gdal.Open(source_file) gt = src_ds.GetGeoTransform() #gt is the spatial resolution #years = range(0,30,1) years = range(0,115,1) #years = range(85,115,1) months = range(1,13,1) yearly_avg = [] for every_year in range(len(years)): yearly_list = 0.0 for month in range(len(months)): for each_coord in range(len(latitude)): index = (years[every_year]*12) + months[month] lat = float(latitude[each_coord]) lon = float(longitude[each_coord]) rb=src_ds.GetRasterBand(index) #rb is the specific raster band plat = int((gt[3] + abs(lat))/gt[1]) plong = int((lon + abs(gt[0]))/gt[1]) if plat > 4477 or plong > 3486: pass else: interval = rb.ReadAsArray(plong,plat,1,1) value = interval[0][0] yearly_list += value/100.0 yearly_avg.append(yearly_list/(len(latitude)*len(months)))#this is where you calculate the average zip_T = list(izip_longest(*[years,yearly_avg])) #with open('CA_Temps_1896_1925.csv','w') as f: #with open('CA_Temps_1996_2015.csv','w') as f: with open('CA_temps_1896_2015.csv','w') as f: writer = csv.writer(f) writer.writerow(['year','temperature']) writer.writerows(zip_T) def sample_through_time(longitude,latitude): 'a function that sample NetCDF files through time for the state of California'\ #(-374495.83635354, 270.0, 0.0, 592636.6658113, 0.0, -270.0) map_ = glob.glob('*.nc') yearly_avg = [] yearly_std = [] yearly_ste = [] Values = [] Months = [] Years = [] Maps = [] Latitude = [] Longitude = [] Sites = [] #years = range(0,30,1) years = range(0,115) #years = range(85,115,1) months = range(1,13,1) for every_year in range(len(years)): yearly_list = [] std_year = [] for maps in range(0,2): source_file = map_[maps] src_ds=gdal.Open(source_file) gt = src_ds.GetGeoTransform() #gt is the spatial resolution for month in range(len(months)): for each_coord in range(len(latitude)): index = (years[every_year]*12) + months[month] lat = float(latitude[each_coord]) lon = float(longitude[each_coord]) rb=src_ds.GetRasterBand(index) #rb is the specific raster band plat = int((gt[3] + abs(lat))/gt[1]) plong = int((lon + abs(gt[0]))/gt[1]) if plat > 4477 or plong > 3486: pass else: interval = rb.ReadAsArray(plong,plat,1,1) value = interval[0][0] if value == -9999: pass else: yearly_list.append(value/100.0) std_year.append(value/100.0) Values.append(value/100.0) Months.append(month+1) Years.append(every_year) Maps.append(maps) Latitude.append(latitude[each_coord]) Longitude.append(longitude[each_coord]) Sites.append(each_coord) yearly_avg.append(numpy.mean(yearly_list)) yearly_std.append(std(std_year)) yearly_ste.append(std(std_year)/sqrt((2.0*len(latitude)*len(months)))) zip_T = list(izip_longest(*[years,yearly_avg,yearly_std,yearly_ste])) #with open('CA_Temps_1896_1925.csv','w') as f: #with open('CA_Temps_1996_2015.csv','w') as f: with open('CA_Temps_1896_2015_TEST.csv','w') as f: writer = csv.writer(f) writer.writerow(['year','temperature','std','ste']) writer.writerows(zip_T) """zip_T = list(izip_longest(*[Values,Months,Years,Maps,Latitude,Longitude,Sites])) with open('single_year_test.csv','w') as f: writer = csv.writer(f) writer.writerow(['Values','Months','Years','Maps','Latitude','Longitude','Sites']) writer.writerows(zip_T)""" def convert_latlon2albers(dataframe): 'convert latitude and longitude to Albers' latitude = make_lat_long_list(dataframe,'latitude') longitude = make_lat_long_list(dataframe,'longitude') epsg3309 = SpatialReference() epsg3309.ImportFromEPSG(3309) epsg4326 = SpatialReference() epsg4326.ImportFromEPSG(4326) latlong2albers = CoordinateTransformation(epsg4326, epsg3309) albers2latlong = CoordinateTransformation(epsg3309,epsg4326) x_conv = [] y_conv = [] for each_coord in range(len(latitude)): new = latlong2albers.TransformPoint(float(longitude[each_coord])*-1,float(latitude[each_coord])) y_conv.append(new[1]) x_conv.append(new[0]) zip_T = list(izip_longest(*[latitude,longitude,x_conv,y_conv])) with open('converted_coordinates.csv','w') as f: writer = csv.writer(f) writer.writerow(['latitude','longitude','latitude_albers','longitude_albers']) writer.writerows(zip_T) #---------------------------# #---------SCRIPT------------# #---------------------------# 'first convert your coordinates to Albers to sample the NetCDF maps' mojave_points = '/path/coordinates.csv' lat_moj = make_lat_long_list(mojave_points,'latitude_albers') lon_moj = make_lat_long_list(mojave_points,'longitude_albers') 'The functions below work to generate a massive dataframe of monthly temperatures for each latitude and longitude' generate_monthly_mins_maxs(directory,lat_moj,lon_moj,full_spp,mass) sample_netCDF_CA_by_year(lat_moj,lon_moj) sample_through_time(lat_moj,lon_moj) 'convert latitude and longitude to Albers CRS' convert_latlon2albers(mojave_points)