forked from NCPP/ocgis
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
2 changed files
with
91 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,46 @@ | ||
rm(list = ls()) | ||
require(plyr) | ||
require(foreign) | ||
|
||
## function for binning | ||
bin = function(dfr,bounds){ | ||
ndays = nrow(dfr) | ||
percents = list() | ||
dfr = subset(dfr,!is.na(VALUE)) | ||
for (name in names(bounds)){ | ||
vec = bounds[[name]] | ||
percents[name] = round((sum(dfr$VALUE >= vec[1] & dfr$VALUE <= vec[2])/ndays)*100,2) | ||
} | ||
return(as.data.frame(percents)) | ||
} | ||
|
||
## some classification variables | ||
bounds = list(caution=c(80,90), | ||
ex_caution=c(91,105), | ||
danger=c(106,130), | ||
ex_danger=c(131,99999)) | ||
|
||
## load time, calculation, user geometry, and value data.frames | ||
tgid = read.csv('tgid.csv') | ||
cid = read.csv('cid.csv') | ||
value = read.csv('value.csv') | ||
|
||
## remove empties and select heat index | ||
value = subset(value,CID == 2) | ||
|
||
## join data to time variables for temporal grouping | ||
value = join(tgid,value,by="TGID") | ||
|
||
## calculate percentages | ||
by_year = ddply(value,c("YEAR","UGID"),bin,bounds) | ||
by_state = ddply(value,c("UGID"),bin,bounds) | ||
|
||
## read the shapefile dbf backing it up first | ||
file.copy('state_boundaries/state_boundaries.dbf','state_boundaries/state_boundaries.dbf.backup') | ||
state_dbf = read.dbf('state_boundaries/state_boundaries.dbf') | ||
|
||
## join the calculated data to it | ||
state_dbf = join(state_dbf,by_state,by="UGID") | ||
|
||
## write it out | ||
write.dbf(state_dbf,'state_boundaries/state_boundaries.dbf') |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,45 @@ | ||
from ocgis.api.interp.iocg.interpreter_ocg import OcgInterpreter | ||
from ocgis.util.shp_cabinet import ShpCabinet | ||
from ocgis.util.inspect import Inspect | ||
|
||
|
||
########## GLOBALS ################################## | ||
## path to nc file | ||
NC = 'TMAXdaysAboveThreshold1950_1955ByGridpoint.nc' | ||
## variable to extract from the nc | ||
VARIABLE = 'gmo_tmax-days_above_threshold' | ||
## directory containing shapefiles | ||
SHP = 'shp' | ||
## name of the folder containing the target shapefile | ||
KEY = 'NC CSC region' | ||
##################################################### | ||
|
||
def main(): | ||
## inspect the dataset | ||
# inspect = Inspect(NC) | ||
|
||
## object to manage shapefile geometries | ||
sc = ShpCabinet(path=SHP) | ||
## extract geometry dictionary from shapefile | ||
geom_dict = sc.get_geom_dict(KEY) | ||
## update id field (required by OCGIS) | ||
for geom in geom_dict: | ||
geom['id'] = geom['CSC_ID'] | ||
## operations dictionary | ||
ops = {'meta':[{'uri':NC,'variable':VARIABLE}], | ||
'time_range':None, | ||
'level_range':None, | ||
'spatial_operation':'intersects', | ||
'output_format':'shp', | ||
'geom':geom_dict, | ||
'aggregate':False} | ||
## execute the operations returning a path to the directory | ||
## containing the output | ||
ret = OcgInterpreter(ops).execute() | ||
|
||
return(ret) | ||
|
||
|
||
if __name__ == '__main__': | ||
ret = main() | ||
print(ret) |