Skip to content

Commit

Permalink
i231: generate demonstration code for geoweb integration
Browse files Browse the repository at this point in the history 
Initial example code for netCDF conversion to GeoJson.
  • Loading branch information
@bekozi
bekozi committed Feb 14, 2014
1 parent 83e115c commit 53b0032
Show file tree
Hide file tree
Showing 2 changed files with 133 additions and 35 deletions.
94 changes: 94 additions & 0 deletions examples/climate_pipes.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
import ocgis
import tempfile
import shutil
import os.path
from ocgis.api.collection import SpatialCollection
from ocgis.conv.fiona_ import GeoJsonConverter

'''
Usage:
- Assumes repo is on current master
- Check the globals below for major script parameters
- The code for generating GeoJson directly from objects is below the code
that using the OcgOperations API. It is commented out.
'''

## path to target file (may also be an OPeNDAP target)
URI = '/home/local/WX/ben.koziol/Dropbox/nesii/project/ocg/bin/climate_data/CanCM4/tas_day_CanCM4_decadal2000_r2i1p1_20010101-20101231.nc'

## the target variable in the dataset to convert
VARIABLE = 'tas'

## this is the path to the shapefile containing state boundaries
SHAPEFILE_PATH = '/home/local/WX/ben.koziol/Dropbox/nesii/project/ocg/bin/shp/state_boundaries/state_boundaries.shp'

## write data to a new temporary directory for each script start
DIR_OUTPUT = tempfile.mkdtemp()

## set the output format to convert to
OUTPUT_FORMAT = 'geojson'
#OUTPUT_FORMAT = 'shp'

## limit the headers in the output.
HEADERS = ['time','year','month','day','value']

## we can either subset the data by a geometry from a shapefile, or convert to
## geojson for the entire spatial domain. there are other options here (i.e. a
## bounding box for tiling or a Shapely geometry).
#GEOM = SHAPEFILE_PATH
GEOM = None


## connect to the dataset and load the data as a field object. this will be used
## to iterate over time coordinates during the conversion step.
rd = ocgis.RequestDataset(uri=URI,variable=VARIABLE)
field = rd.get()

## selecting specific geometries from a shapefile requires knowing the target
## geometry's UGID inside the shapefile. shapefile are required to have this
## identifier at this time as a full attribute search is not enabled. this code
## searches for TX to find the UGID associated with that state.
if GEOM is None:
select_ugid = None
else:
select_geom = filter(lambda x: x['properties']['STATE_ABBR'] == 'TX',ocgis.ShpCabinetIterator(path=SHAPEFILE_PATH))
## this argument must always come in as a list
select_ugid = [select_geom[0]['properties']['UGID']]

## USE OCGIS OPERATIONS ########################################################

## get an example time coordinate
centroid = field.temporal.value_datetime[0]
print('writing geojson for time slice: {0}'.format(centroid))
## this again is the target dataset with a time range for subsetting now
rd = ocgis.RequestDataset(uri=URI,variable=VARIABLE,time_range=[centroid,centroid])
## name of the output geojson file
prefix = 'ocgis_output_{0}'.format(OUTPUT_FORMAT)
## parameterize the operations to be performed on the target dataset
ops = ocgis.OcgOperations(dataset=rd,geom=SHAPEFILE_PATH,select_ugid=select_ugid,
output_format=OUTPUT_FORMAT,prefix=prefix,dir_output=DIR_OUTPUT,
headers=HEADERS)
ret = ops.execute()
print('path to output file: {0}'.format(ret))

## USE OCGIS OBJECTS DIRECTLY ##################################################

### get an example time coordinate
#centroid = field.temporal.value_datetime[0]
#print('writing geojson for time slice: {0}'.format(centroid))
#
### subset the field for that coordinate
#subsetted = field.get_between('temporal',lower=centroid,upper=centroid)
#
### alternatively you may simple slice the data (5-d slicing only - realization,
### time,level,row,column
##subsetted = field[:,0,:,:,:]
#
### collection objects are used the converters
#sc = SpatialCollection(headers=HEADERS)
### add the field. a geometry identifier value is required (GIS first always!).
### the geometry is empty however. the field needs to be given an alias for unique
### storage in the collection. the last argument is field object itself.
#sc.add_field(1,None,'tas',subsetted)
#conv = GeoJsonConverter([sc],DIR_OUTPUT,'ocgis_output_geojson')
#print(conv.write())
74 changes: 39 additions & 35 deletions src/ocgis/conv/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ def write(self):

try:
build = True
if self._add_ugeom and self.ops.geom is not None:
if self._add_ugeom and self.ops is not None and self.ops.geom is not None:
write_geom = True
else:
write_geom = False
Expand Down Expand Up @@ -155,42 +155,46 @@ def write(self):
except UnboundLocalError:
pass

## added OCGIS metadata output if requested.
if self.add_meta:
ocgis_lh('adding OCGIS metadata file','conv',logging.DEBUG)
lines = MetaConverter(self.ops).write()
out_path = os.path.join(self.outdir,self.prefix+'_'+MetaConverter._meta_filename)
with open(out_path,'w') as f:
f.write(lines)

## add the dataset descriptor file if specified
if self._add_did_file:
ocgis_lh('writing dataset description (DID) file','conv',logging.DEBUG)
from ocgis.conv.csv_ import OcgDialect
## the metadata and dataset descriptor files may only be written if
## OCGIS operations are present.
if self.ops is not None:
## added OCGIS metadata output if requested.
if self.add_meta:
ocgis_lh('adding OCGIS metadata file','conv',logging.DEBUG)
lines = MetaConverter(self.ops).write()
out_path = os.path.join(self.outdir,self.prefix+'_'+MetaConverter._meta_filename)
with open(out_path,'w') as f:
f.write(lines)

headers = ['DID','VARIABLE','ALIAS','URI','STANDARD_NAME','UNITS','LONG_NAME']
out_path = os.path.join(self.outdir,self.prefix+'_did.csv')
with open(out_path,'w') as f:
writer = csv.writer(f,dialect=OcgDialect)
writer.writerow(headers)
## add the dataset descriptor file if specified and OCGIS operations
## are present.
if self._add_did_file:
ocgis_lh('writing dataset description (DID) file','conv',logging.DEBUG)
from ocgis.conv.csv_ import OcgDialect

headers = ['DID','VARIABLE','ALIAS','URI','STANDARD_NAME','UNITS','LONG_NAME']
out_path = os.path.join(self.outdir,self.prefix+'_did.csv')
with open(out_path,'w') as f:
writer = csv.writer(f,dialect=OcgDialect)
writer.writerow(headers)
for rd in self.ops.dataset:
row = [rd.did,rd.variable,rd.alias,rd.uri]
ref_variable = rd._source_metadata['variables'][rd.variable]['attrs']
row.append(ref_variable.get('standard_name',None))
row.append(ref_variable.get('units',None))
row.append(ref_variable.get('long_name',None))
writer.writerow(row)

## add source metadata if requested
if self._add_source_meta:
ocgis_lh('writing source metadata file','conv',logging.DEBUG)
out_path = os.path.join(self.outdir,self.prefix+'_source_metadata.txt')
to_write = []
for rd in self.ops.dataset:
row = [rd.did,rd.variable,rd.alias,rd.uri]
ref_variable = rd._source_metadata['variables'][rd.variable]['attrs']
row.append(ref_variable.get('standard_name',None))
row.append(ref_variable.get('units',None))
row.append(ref_variable.get('long_name',None))
writer.writerow(row)

## add source metadata if requested
if self._add_source_meta:
ocgis_lh('writing source metadata file','conv',logging.DEBUG)
out_path = os.path.join(self.outdir,self.prefix+'_source_metadata.txt')
to_write = []
for rd in self.ops.dataset:
ip = Inspect(request_dataset=rd)
to_write += ip.get_report()
with open(out_path,'w') as f:
f.writelines('\n'.join(to_write))
ip = Inspect(request_dataset=rd)
to_write += ip.get_report()
with open(out_path,'w') as f:
f.writelines('\n'.join(to_write))

## return the internal path unless overloaded by subclasses.
ret = self._get_return_()
Expand Down

0 comments on commit 53b0032

Please sign in to comment.