diff --git a/combine_features.py b/combine_features.py index b317e571..8c998a42 100755 --- a/combine_features.py +++ b/combine_features.py @@ -1,14 +1,14 @@ #!/usr/bin/env python """ -This script takes a file containing one or more feature definitions, pointed +This script takes a file containing one or more feature definitions, pointed to by the -f flag, and a new name for the combined feature, provided through the -n flag. The geometry of the feature definitions are combined into a single feature definition, which is placed in (or appended to) the file pointed to with the -o flag (features.geojson by default). Author: Xylar Asay-Davis -Last Modified: 09/29/2016 +Last Modified: 10/16/2016 """ import json @@ -32,9 +32,6 @@ parser.add_argument("-n", "--new_feature_name", dest="new_feature_name", help="The new name of the combined feature", metavar="NAME", required=True) -parser.add_argument("-g", "--groupName", dest="groupName", - help="Feature group name.", - metavar="GROUPNAME", default="unspecifiedGroupName") parser.add_argument("-o", "--output", dest="output_file_name", help="Output file, e.g., features.geojson.", metavar="PATH", default="features.geojson") @@ -100,7 +97,6 @@ feature['properties']['constituents'] = '; '.join(list(set(featureNames))) feature['geometry'] = shapely.geometry.mapping(combinedShape) features['features'].append(feature) -features['groupName'] = args.groupName if feature['geometry']['type'] == 'GeometryCollection': print "Error: combined geometry from %s is of type GeometryCollection."%(args.feature_file) diff --git a/difference_features.py b/difference_features.py index b661d43a..73f1f5e9 100755 --- a/difference_features.py +++ b/difference_features.py @@ -1,14 +1,14 @@ #!/usr/bin/env python """ This script takes a file containing one or more feature definitions, that is -pointed to by the -f flag and a second set of one or more masking feature +pointed to by the -f flag and a second set of one or more masking feature definition, pointed to with the -m flag. The masking features are masked out -of (i.e. removed from) the original feature definitions. The resulting -features are placed in (or appended to) the output file pointed to with the +of (i.e. removed from) the original feature definitions. The resulting +features are placed in (or appended to) the output file pointed to with the -o flag (features.geojson by default). Authors: Xylar Asay-Davis -Last Modified: 09/29/2016 +Last Modified: 10/16/2016 """ import json @@ -31,9 +31,6 @@ help="Feature file with one or more features whose overlap " "with features in feature_file should be removed", metavar="FILE2", required=True) -parser.add_argument("-g", "--groupName", dest="groupName", - help="Feature group name.", - metavar="GROUPNAME", default="unspecifiedGroupName") parser.add_argument("-o", "--output", dest="output_file_name", help="Output file, e.g., features.geojson.", metavar="PATH", default="features.geojson") @@ -111,8 +108,6 @@ else: print "%s has been removed."%name -features['groupName'] = args.groupName - write_all_features(features, out_file_name, indent=4) # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/driver_scripts/setup_MOC_basins.py b/driver_scripts/setup_MOC_basins.py deleted file mode 100755 index d69e4203..00000000 --- a/driver_scripts/setup_MOC_basins.py +++ /dev/null @@ -1,67 +0,0 @@ -#!/usr/bin/env python -""" -This script creates basins for computing the meridional overturning -circulation (MOC). No arguments are required. The optional --plot -flag can be used to produce plots of each MOC basin. - -Authors: Xylar Asay-Davis, Phillip J. Wolfram -Last Modified: 09/09/2016 -""" -import os -import subprocess -from optparse import OptionParser - -def spcall(args): #{{{ - return subprocess.check_call(args, env=os.environ.copy()) #}}} - -parser = OptionParser() -parser.add_option("--plot", action="store_true", dest="plot") - -options, args = parser.parse_args() - -groupNameMOC = 'MOCBasinRegionsGroup' -groupNameBasins = 'OceanBasinRegionsMOCGroup' - -subBasins = {'Atlantic': ['Atlantic', 'Mediterranean'], - 'IndoPacific': ['Pacific', 'Indian'], - 'Pacific': ['Pacific'], - 'Indian': ['Indian']} - -MOCMaskFileNames = {'Atlantic': 'ocean/region/MOC_mask_34S/region.geojson', - 'IndoPacific': 'ocean/region/MOC_mask_34S/region.geojson', - 'Pacific': 'ocean/region/MOC_mask_6S/region.geojson', - 'Indian': 'ocean/region/MOC_mask_6S/region.geojson'} - -for basinName in subBasins: - - # output file names - basinFileName = '%s_Basin_separate.geojson'%basinName - basinCombinedFileName = '%s_Basin.geojson'%basinName - MOCName = '%s_MOC.geojson'%basinName - imageName = '%s_MOC.png'%basinName - - # remove old files (to ensure there isn't double-appending via merge_features - for afile in [basinFileName, basinCombinedFileName, MOCName, imageName]: - if os.path.exists(afile): - os.remove(afile) - - print " * merging features to make %s Basin"%basinName - - for oceanName in subBasins[basinName]: - spcall(['./merge_features.py', '-d', 'ocean/region', '-t', '%s_Basin'%oceanName, - '-o', basinFileName]) - - #merge the the features into a single file - print " * combining features into single feature named %s_MOC"%basinName - spcall(['./combine_features.py', '-f', basinFileName, '-n', '%s_MOC'%basinName, - '-g', groupNameBasins, '-o', basinCombinedFileName]) - - print " * masking out features south of MOC region" - spcall(['./difference_features.py', '-f', basinCombinedFileName, - '-m', MOCMaskFileNames[basinName], '-g', groupNameMOC, '-o', MOCName]) - - if options.plot: - spcall(['./plot_features.py', '-f', MOCName, '-o', imageName, '-m', 'cyl']) - - -# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/driver_scripts/setup_ocean_basins.py b/driver_scripts/setup_ocean_basins.py index 040f7943..7c7fd94c 100755 --- a/driver_scripts/setup_ocean_basins.py +++ b/driver_scripts/setup_ocean_basins.py @@ -1,6 +1,6 @@ #!/usr/bin/env python """ -This script creates 6 major ocean basins. No arguments are required. The +This script creates 6 major ocean basins. No arguments are required. The optional --plot flag can be used to produce plots of each basin. """ import os @@ -12,35 +12,165 @@ def spcall(args): #{{{ return subprocess.check_call(args, env=os.environ.copy()) #}}} +def build_basins(): #{{{ + basinGroupName = 'OceanBasinRegionsGroup' + basinFileName = 'oceanBasins.geojson' + + # temp file names that we delete later + tempSeparateBasinFileName = 'temp_separate_basins.geojson' + tempCombinedBasinFileName = 'temp_combined_basins.geojson' + + # remove old files so we don't unintentionally append features + for fileName in [basinFileName, + tempSeparateBasinFileName, + tempCombinedBasinFileName]: + if os.path.exists(fileName): + os.remove(fileName) + + # build ocean basins from regions with the appropriate tags + for oceanName in ['Atlantic', 'Pacific', 'Indian', 'Arctic', 'Southern_Ocean', + 'Mediterranean']: + + tag = '%s_Basin'%oceanName + + print " * merging features to make %s Basin"%oceanName + spcall(['./merge_features.py', '-d', 'ocean/region', '-t', tag, + '-o', tempSeparateBasinFileName]) + + #merge the the features into a single file + print " * combining features into single feature named %s_Basin"%oceanName + spcall(['./combine_features.py', '-f', tempSeparateBasinFileName, + '-n', '%s_Basin'%oceanName, + '-o', tempCombinedBasinFileName]) + + if(options.plot): + spcall(['./plot_features.py', '-f', tempCombinedBasinFileName, + '-o', '%s_Basin.png'%oceanName, '-m', 'cyl']) + + spcall(['./merge_features.py', '-f', tempCombinedBasinFileName, + '-o', basinFileName]) + + # remove temp files + for fileName in [tempSeparateBasinFileName, tempCombinedBasinFileName]: + os.remove(fileName) + + + spcall(['./set_group_name.py', '-f', basinFileName, + '-g', basinGroupName]) + + if(options.plot): + spcall(['./plot_features.py', '-f', basinFileName, + '-o', 'oceanBasins.png', '-m', 'cyl']) #}}} + + +def build_MOC_basins(): #{{{ + MOCFileName = 'MOCBasins.geojson' + MOCGroupName = 'MOCBasinRegionsGroup' + + MOCSouthernBoundaryFileName = 'MOCSounthernBoundaries.geojson' + MOCSouthernBoundaryGroupName = 'MOCSouthernBoundaryGroup' + + MOCSubBasins = {'Atlantic': ['Atlantic', 'Mediterranean'], + 'IndoPacific': ['Pacific', 'Indian'], + 'Pacific': ['Pacific'], + 'Indian': ['Indian']} + + + MOCSouthernBoundary = {'Atlantic': '34S', + 'IndoPacific': '34S', + 'Pacific': '6S', + 'Indian': '6S'} + + # temp file names that we delete later + tempSeparateBasinFileName = 'temp_separate_basins.geojson' + tempCombinedBasinFileName = 'temp_combined_basins.geojson' + tempMOCFileName = 'temp_MOC.geojson' + tempSouthernBoundaryFileName = 'temp_southern_boundary.geojson' + tempSouthernBoundaryRenamedFileName = 'temp_southern_boundary_renamed.geojson' + + # remove old files so we don't unintentionally append features + for fileName in [MOCFileName, + MOCSouthernBoundaryFileName, + tempSeparateBasinFileName, + tempCombinedBasinFileName, + tempMOCFileName, + tempSouthernBoundaryFileName, + tempSouthernBoundaryRenamedFileName]: + if os.path.exists(fileName): + os.remove(fileName) + + + # build MOC basins + for basinName in MOCSubBasins: + + imageName = '%s_MOC.png'%basinName + + MOCMaskFileName = 'ocean/region/MOC_mask_%s/region.geojson' \ + % MOCSouthernBoundary[basinName] + + MOCSouternBoundaryFileName = 'ocean/transect/MOC_%s/transect.geojson' \ + % MOCSouthernBoundary[basinName] + + print " * merging features to make %s Basin"%basinName + + for oceanName in MOCSubBasins[basinName]: + spcall(['./merge_features.py', '-d', 'ocean/region', + '-t', '%s_Basin'%oceanName, + '-o', tempSeparateBasinFileName]) + + #merge the the features into a single file + print " * combining features into single feature named %s_MOC"%basinName + spcall(['./combine_features.py', '-f', tempSeparateBasinFileName, + '-n', '%s_MOC'%basinName, + '-o', tempCombinedBasinFileName]) + + print " * masking out features south of MOC region" + spcall(['./difference_features.py', '-f', tempCombinedBasinFileName, + '-m', MOCMaskFileName, + '-o', tempMOCFileName]) + + spcall(['./merge_features.py', '-f', tempMOCFileName, + '-o', MOCFileName]) + + if options.plot: + spcall(['./plot_features.py', '-f', tempMOCFileName, '-o', imageName, + '-m', 'cyl']) + + print " * make feature for southern boundary of MOC region" + spcall(['./intersect_features.py', '-f', tempCombinedBasinFileName, + '-i', MOCSouternBoundaryFileName, + '-o', tempSouthernBoundaryFileName]) + + # use combine_features to rename the feature + spcall(['./combine_features.py', '-f', tempSouthernBoundaryFileName, + '-n', '%s_MOC_southern_boundary'%basinName, + '-o', tempSouthernBoundaryRenamedFileName]) + + spcall(['./merge_features.py', '-f', tempSouthernBoundaryRenamedFileName, + '-o', MOCSouthernBoundaryFileName]) + + # remove temp files + for fileName in [tempSeparateBasinFileName, tempCombinedBasinFileName, + tempMOCFileName, tempSouthernBoundaryFileName, + tempSouthernBoundaryRenamedFileName]: + os.remove(fileName) + + spcall(['./set_group_name.py', '-f', MOCFileName, + '-g', MOCGroupName]) + + spcall(['./set_group_name.py', '-f', MOCSouthernBoundaryFileName, + '-g', MOCSouthernBoundaryGroupName]) + + if options.plot: + spcall(['./plot_features.py', '-f', MOCFileName, + '-o', 'MOCBasins.png', '-m', 'cyl']) parser = OptionParser() parser.add_option("--plot", action="store_true", dest="plot") options, args = parser.parse_args() -groupName = 'OceanBasinRegionsGroup' - -for oceanName in 'Atlantic', 'Pacific', 'Indian', 'Arctic', 'Southern_Ocean', 'Mediterranean': - - tag = '%s_Basin'%oceanName - basinFileNameName = '%s_Basin_separate.geojson'%oceanName - basinCombinedFileName = '%s_Basin.geojson'%oceanName - imageName = '%s_Basin.png'%oceanName - - # remove old files (to ensure there isn't double-appending via merge_features - for afile in [basinFileNameName, basinCombinedFileName, imageName]: - if os.path.exists(afile): - os.remove(afile) - - print " * merging features to make %s Basin"%oceanName - args = ['./merge_features.py', '-d', 'ocean/region', '-t', tag, '-o', basinFileNameName] - subprocess.check_call(args, env=os.environ.copy()) - - #merge the the features into a single file - print " * combining features into single feature named %s_Basin"%oceanName - spcall(['./combine_features.py', '-f', basinFileNameName, '-n', '%s_Basin'%oceanName, - '-g', groupName, '-o', basinCombinedFileName]) - - if(options.plot): - args = ['./plot_features.py', '-f', basinCombinedFileName, '-o', imageName, '-m', 'cyl'] - subprocess.check_call(args, env=os.environ.copy()) +build_basins() +build_MOC_basins() + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/intersect_features.py b/intersect_features.py new file mode 100755 index 00000000..254e412f --- /dev/null +++ b/intersect_features.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python +""" +This script takes a file containing one or more feature definitions, +indicated by the -f flag, and a second set of one or more intersecting feature +definition, pointed to with the -i flag. The resulting feature definitions +are the intersection of the original feature definitions with the all of the +intersecting features (i.e. the intersection of the intersecting features). +The resulting features are placed in (or appended to) the output file indicated +with the -o flag (features.geojson by default). + +Authors: Xylar Asay-Davis +Last Modified: 10/16/2016 +""" + +import json +import argparse +from collections import defaultdict +from utils.feature_write_utils import write_all_features +from utils.feature_test_utils import feature_already_exists + +import os.path + +import shapely.geometry +import shapely.ops + +parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawTextHelpFormatter) +parser.add_argument("-f", "--feature_file", dest="feature_file", + help="Feature file to be intersected", metavar="FILE1", + required=True) +parser.add_argument("-i", "--intersection_file", dest="intersection_file", + help="Feature file with features whose intersection with" + "those in feature_file should be retained", + metavar="FILE2", required=True) +parser.add_argument("-g", "--groupName", dest="groupName", + help="Group name for output features", + metavar="GROUPNAME", default="enterGroupName") +parser.add_argument("-o", "--output", dest="output_file_name", + help="Output file name", + metavar="PATH", default="features.geojson") + +args = parser.parse_args() + + +out_file_name = args.output_file_name + +features = defaultdict(list) + +if os.path.exists(out_file_name): + try: + with open(out_file_name) as f: + appended_file = json.load(f) + for feature in appended_file['features']: + features['features'].append(feature) + del appended_file + except: + pass + +featuresToIntersect = defaultdict(list) + +try: + with open(args.feature_file) as f: + feature_file = json.load(f) + + for feature in feature_file['features']: + featuresToIntersect['features'].append(feature) + + del feature_file +except: + print "Error parsing geojson file: %s" % (args.feature_file) + raise + +intersectionShape = None +try: + with open(args.intersection_file) as f: + feature_file = json.load(f) + + for feature in feature_file['features']: + shape = shapely.geometry.shape(feature['geometry']) + if intersectionShape is None: + intersectionShape = shape + else: + intersectionShape = shape.intersection(intersectionShape) + + del feature_file +except: + print "Error parsing geojson file: %s" % (args.mask_file) + raise + +for feature in featuresToIntersect['features']: + name = feature['properties']['name'] + if feature_already_exists(features, feature): + print "Warning: feature %s already in features.geojson. " \ + "Skipping..." % name + continue + featureShape = shapely.geometry.shape(feature['geometry']) + if featureShape.intersects(intersectionShape): + print "%s has been intersected." % name + featureShape = featureShape.intersection(intersectionShape) + assert(not featureShape.is_empty) + feature['geometry'] = shapely.geometry.mapping(featureShape) + features['features'].append(feature) + else: + print "Warning: feature %s has no intersection with features in %s." \ + "Skipping..." % (name, args.intersection_file) + +features['groupName'] = args.groupName + +write_all_features(features, out_file_name, indent=4) + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/merge_features.py b/merge_features.py index 820d1c95..391551fe 100755 --- a/merge_features.py +++ b/merge_features.py @@ -17,12 +17,13 @@ input line. Authors: Douglas Jacobsen, Xylar Asay-Davis -Last Modified: 9/29/2016 +Last Modified: 10/16/2016 """ -import sys, os, glob, shutil, numpy, fnmatch +import os import json import argparse +import fnmatch from collections import defaultdict from utils.feature_write_utils import write_all_features from utils.feature_test_utils import match_tag_list, feature_already_exists @@ -30,7 +31,6 @@ parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawTextHelpFormatter) parser.add_argument("-f", "--feature_file", dest="feature_file", help="Single feature file to append to features.geojson", metavar="FILE") parser.add_argument("-d", "--features_directory", dest="features_dir", help="Directory containing multiple feature files, each will be appended to features.geojson", metavar="PATH") -parser.add_argument("-g", "--group", dest="groupName", help="Feature group name.", metavar="GROUPNAME", default="unspecifiedGroupName") parser.add_argument("-t", "--tags", dest="tags", help="Semicolon separated list of tags to match features against.", metavar='"TAG1;TAG2;TAG3"') parser.add_argument("-o", "--output", dest="output_file_name", help="Output file, e.g., features.geojson.", metavar="PATH", default="features.geojson") @@ -102,8 +102,6 @@ print "Error parsing geojson file: %s"%(path) del paths -all_features['groupName'] = args.groupName - write_all_features(all_features, file_to_append, indent=4) # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/ocean/transect/MOC_34S/transect.geojson b/ocean/transect/MOC_34S/transect.geojson new file mode 100644 index 00000000..efa094af --- /dev/null +++ b/ocean/transect/MOC_34S/transect.geojson @@ -0,0 +1,32 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "object": "transect", + "component": "ocean", + "name": "MOC 34S", + "tags": "MOC_34S", + "author": "Xylar Asay-Davis" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + -180.000000, + -34.000000 + ], + [ + 0.000000, + -34.000000 + ], + [ + 180.000000, + -34.000000 + ] + ] + } + } + ] +} \ No newline at end of file diff --git a/ocean/transect/MOC_6S/transect.geojson b/ocean/transect/MOC_6S/transect.geojson new file mode 100644 index 00000000..c0b43420 --- /dev/null +++ b/ocean/transect/MOC_6S/transect.geojson @@ -0,0 +1,32 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "object": "transect", + "component": "ocean", + "name": "MOC 6S", + "tags": "MOC_6S", + "author": "Xylar Asay-Davis" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + -180.000000, + -6.000000 + ], + [ + 0.000000, + -6.000000 + ], + [ + 180.000000, + -6.000000 + ] + ] + } + } + ] +} \ No newline at end of file diff --git a/set_group_name.py b/set_group_name.py new file mode 100755 index 00000000..d70c2640 --- /dev/null +++ b/set_group_name.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python +""" +This script reads features from the input features file pointed to with the -f +flag, adds a group name given with the -g flag, and writes the features to the +same file. + +Authors: Xylar Asay-Davis +Last Modified: 10/16/2016 +""" + +import os +import json +import argparse +from collections import defaultdict +from utils.feature_write_utils import write_all_features +from utils.feature_test_utils import feature_already_exists + +parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawTextHelpFormatter) +parser.add_argument("-f", "--feature_file", dest="feature_file", required=True, + help="Input and output feature file where group name is to" + "be set", metavar="FILE") +parser.add_argument("-g", "--group", dest="groupName", + help="Feature group name", metavar="GROUPNAME", + required=True) + +args = parser.parse_args() + +if not os.path.exists(args.feature_file): + parser.error('The file %s does not exist.'%(args.feature_file)) + +all_features = defaultdict(list) + +with open(args.feature_file) as f: + feature_file = json.load(f) + + for feature in feature_file['features']: + if not feature_already_exists(all_features, feature): + all_features['features'].append(feature) + +all_features['groupName'] = args.groupName + +write_all_features(all_features, args.feature_file, indent=4) + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/split_features.py b/split_features.py index addc0684..8e50043b 100755 --- a/split_features.py +++ b/split_features.py @@ -10,7 +10,7 @@ to determine the base directory. Authors: Douglas Jacobsen, Xylar Asay-Davis -Last Modified: 9/29/2016 +Last Modified: 10/16/2016 """ @@ -57,6 +57,7 @@ out_file_name = '%s/%s/%s/%s.geojson'%(base_dir, object_type, dir_name, object_type) - write_all_features({'features':[feature]}, out_file_name, indent=4) + write_all_features({'features':[feature]}, out_file_name, indent=4, + defaultGroupName=None) # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/utils/feature_write_utils.py b/utils/feature_write_utils.py index 44418d9f..66288858 100644 --- a/utils/feature_write_utils.py +++ b/utils/feature_write_utils.py @@ -13,7 +13,8 @@ from collections import OrderedDict -def write_all_features(features, out_file_name, indent=4): # {{{ +def write_all_features(features, out_file_name, indent=4, + defaultGroupName='enterGroupName'): # {{{ json.encoder.FLOAT_REPR = lambda o: format(o, 'f') for index in range(len(features['features'])): @@ -27,6 +28,8 @@ def write_all_features(features, out_file_name, indent=4): # {{{ outFeatures = OrderedDict((('type', features['type']),)) if 'groupName' in features.keys(): outFeatures['groupName'] = features['groupName'] + elif defaultGroupName is not None: + outFeatures['groupName'] = defaultGroupName # Add the rest for key in features: