updated v0.2.7

- Improved overlap calculations for larger geometries - Added a geometry simplification tool to reduce number of vertices
tyson-swetnam · Apr 25, 2019 · f3ec557 · f3ec557
1 parent 313b4fd
commit f3ec557
Show file tree

Hide file tree

Showing 9 changed files with 136 additions and 21 deletions.
diff --git a/README.md b/README.md
@@ -21,6 +21,7 @@ http://doi.org/10.5281/zenodo.2635407
 * [porder Ordersv2 Simple Client](#porder-ordersv2-simple-client)
     * [porder quota](#porder-quota)
     * [base64](#base64)
+    * [simplify](#simplify)
     * [idlist](#idlist)
     * [difflist](#difflist)
     * [idsplit](#idsplit)
@@ -42,7 +43,7 @@ Shapely is notoriously difficult as a library to install on windows machines so
 
 ```pip install Shapely-1.6.4.post1-cp27-cp27m-win32.whl```
 
-Or you can use [anaconda to install](https://conda-forge.github.io/). Again, both of these options are mentioned on [Shapely’s Official PyPI page](https://pypi.org/project/Shapely/)
+Or you can use [anaconda to install](https://conda-forge.github.io/). Again, both of these options are mentioned on [Shapely’s Official PyPI page](https://pypi.org/project/Shapely/). **Fiona** is a recommended install used by the simplify tool, but it is not necessary. You can find installation instructions [here](https://pypi.org/project/Fiona/1.8.6/#description)
 
 Once you have shapely configured. To install **porder: Simple CLI for Planet ordersv2 API** you can install using two methods
 
@@ -68,7 +69,7 @@ Installation is an optional step; the application can be also run directly by ex
 
 Make sure you initialized planet client by typing ```planet init``` or ```export``` or ```set PL_API_KEY=Your API Key``` As usual, to print help:
 
-![porder_cli](https://user-images.githubusercontent.com/28806922/53095363-a1059500-34ea-11e9-840f-8e155e3d0ade.png)
+![porder_cli](https://user-images.githubusercontent.com/6677629/56763897-75fd2300-6771-11e9-97b7-97a4780c81c9.png)
 
 To obtain help for a specific functionality, simply call it with _help_ switch, e.g.: `porder idlist -h`. If you didn't install porder, then you can run it just by going to *porder* directory and running `python porder.py [arguments go here]`
 
@@ -85,6 +86,11 @@ This does exactly as it sounds, it encodes your credential files to base64 for u
 
 ![porder_base64](https://user-images.githubusercontent.com/28806922/53096754-d495ee80-34ed-11e9-980d-418601bc975a.png)
 
+### simplify
+This reduces the number of vertices for a multi vertex and complex GeoJSON. Extremely high vertex count (over 500) seem to fail and hence this tool allows you to export a new geojson with reduced vertices. It uses an implementation of the Visvalingam-Wyatt line simplification algorithm. This tool does work with and without Fiona, but Fiona installation is recommended.
+
+![porder simplify](https://user-images.githubusercontent.com/6677629/56763793-36ced200-6771-11e9-8b61-8f94b1f61152.png)
+
 ### idlist
 Create an idlist for your geometry based on some basic filters,including geometry, start and end date and cloud cover. If no cloud cover is specified everything form 0 to 100% cloud cover is included. For now the tool can handle geojson,json and kml files. The output is a csv file with ids. The tool also allows you to make sure you get percentage overlap, when selecting image, for clip operations adjust it accordingly (usally --ovp 1 for orders not to fail during clip). The tool now also prints estimated area in Square kilometes for the download and estimated area if you clipped your area with the geometry you are searching (just estimates).
 
@@ -193,6 +199,10 @@ A simple setup would be
 
 ## Changelog
 
+### v0.2.7
+- Improved overlap calculations for larger geometries
+- Added a geometry simplification tool to reduce number of vertices
+
 ### v0.2.6
 - Skysat area are calculated using EPSG:3857 to resolve metadata EPSG issue
 - General improvements

diff --git a/dist/porder-0.2.6.tar.gz b/dist/porder-0.2.6.tar.gz
diff --git a/dist/porder-0.2.7.tar.gz b/dist/porder-0.2.7.tar.gz
diff --git a/porder/__init__.py b/porder/__init__.py
@@ -2,4 +2,4 @@
 
 __author__ = 'Samapriya Roy'
 __email__ = '[email protected]'
-__version__ = '0.2.6'
+__version__ = '0.2.7'
diff --git a/porder/bundles.json b/porder/bundles.json
@@ -22,6 +22,7 @@
         "PSOrthoTile": [
           "analytic",
           "udm",
+          "udm2",
           "analytic_xml"
         ],
         "PSScene3Band": [

diff --git a/porder/geojson2id.py b/porder/geojson2id.py
@@ -48,6 +48,15 @@
 far=[]
 n=0
 
+#get coordinates list depth
+def list_depth(dic, level = 1):
+    counter=0
+    str_dic = str(dic)
+    if "[[[[" in str_dic:
+        counter += 1
+    return(counter)
+
+
 def handle_page(page,item,asset,num,outfile,gmain,ovp):
     global n
     if num is None:
@@ -73,13 +82,20 @@ def handle_page(page,item,asset,num,outfile,gmain,ovp):
                             #print('ID '+str(it)+' has percentage overlap: '+str(intersect.area/geommain.area*100))
                             proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                                 pyproj.Proj(init='epsg:'+str(epsgcode)))
-                            if (intersect.area/gmain.area)*100>=ovp:
-                                ar.append(transform(proj,intersect).area/1000000)
-                                far.append(transform(proj,geom2).area/1000000)
-                                # print('ID '+str(it)+' has percentage overlap: '+str(intersect.area/geom2.area*100))
-                                with open(outfile,'a') as csvfile:
-                                    writer=csv.writer(csvfile,delimiter=',',lineterminator='\n')
-                                    writer.writerow([it])
+                            if transform(proj,gmain).area>transform(proj,geom2).area:
+                                if (transform(proj,intersect).area/transform(proj,geom2).area*100)>=ovp:
+                                    ar.append(transform(proj,intersect).area/1000000)
+                                    far.append(transform(proj,geom2).area/1000000)
+                                    with open(outfile,'a') as csvfile:
+                                        writer=csv.writer(csvfile,delimiter=',',lineterminator='\n')
+                                        writer.writerow([it])
+                            elif transform(proj,geom2).area>transform(proj,gmain).area:
+                                if (transform(proj,intersect).area/transform(proj,gmain).area*100)>=ovp:
+                                    ar.append(transform(proj,intersect).area/1000000)
+                                    far.append(transform(proj,geom2).area/1000000)
+                                    with open(outfile,'a') as csvfile:
+                                        writer=csv.writer(csvfile,delimiter=',',lineterminator='\n')
+                                        writer.writerow([it])
         except Exception as e:
             print(e)
     elif num is not None:
@@ -106,15 +122,26 @@ def handle_page(page,item,asset,num,outfile,gmain,ovp):
                             #print('ID '+str(it)+' has percentage overlap: '+str(intersect.area/geommain.area*100))
                             proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                                 pyproj.Proj(init='epsg:'+str(epsgcode)))
-                            if (intersect.area/gmain.area)*100>=ovp:
-                                ar.append(transform(proj,intersect).area/1000000)
-                                far.append(transform(proj,geom2).area/1000000)
-                                # print('ID '+str(it)+' has percentage overlap: '+str(intersect.area/geom2.area*100))
-                                n=n+1
-                                #print(n)
-                                with open(outfile,'a') as csvfile:
-                                    writer=csv.writer(csvfile,delimiter=',',lineterminator='\n')
-                                    writer.writerow([it])
+                            if transform(proj,gmain).area>transform(proj,geom2).area:
+                                if (transform(proj,intersect).area/transform(proj,geom2).area*100)>=ovp:
+                                    ar.append(transform(proj,intersect).area/1000000)
+                                    far.append(transform(proj,geom2).area/1000000)
+                                    # print('ID '+str(it)+' has percentage overlap: '+str(intersect.area/geom2.area*100))
+                                    n=n+1
+                                    #print(n)
+                                    with open(outfile,'a') as csvfile:
+                                        writer=csv.writer(csvfile,delimiter=',',lineterminator='\n')
+                                        writer.writerow([it])
+                            elif transform(proj,geom2).area>transform(proj,gmain).area:
+                                if (transform(proj,intersect).area/transform(proj,gmain).area*100)>=ovp:
+                                    ar.append(transform(proj,intersect).area/1000000)
+                                    far.append(transform(proj,geom2).area/1000000)
+                                    # print('ID '+str(it)+' has percentage overlap: '+str(intersect.area/geom2.area*100))
+                                    n=n+1
+                                    #print(n)
+                                    with open(outfile,'a') as csvfile:
+                                        writer=csv.writer(csvfile,delimiter=',',lineterminator='\n')
+                                        writer.writerow([it])
             data=csv.reader(open(outfile).readlines()[0: num])
 
             with open(outfile, "w") as f:
@@ -149,7 +176,12 @@ def idl(infile,start,end,item,asset,num,cmin,cmax,outfile,ovp):
         if infile.endswith('.geojson'):
             with open(infile) as aoi:
                 aoi_resp = json.load(aoi)
-                aoi_geom = aoi_resp['features'][0]['geometry']['coordinates']
+                if list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==0:
+                    aoi_geom = aoi_resp['features'][0]['geometry']['coordinates']
+                elif list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==1:
+                    aoi_geom = aoi_resp['features'][0]['geometry']['coordinates'][0]
+                else:
+                    print('Please check GeoJSON: Could not parse coordinates')
         elif infile.endswith('.json'):
             with open (infile) as aoi:
                 aoi_resp=json.load(aoi)

diff --git a/porder/geojson_simplify.py b/porder/geojson_simplify.py
@@ -0,0 +1,60 @@
+import visvalingamwyatt as vw
+import json
+
+
+# get coordinates list depth
+def list_depth(dic, level=1):
+    counter = 0
+    str_dic = str(dic)
+    if "[[[[" in str_dic:
+        counter += 1
+    return(counter)
+
+
+def geosimple(inp,output,num):
+    try:
+        import fiona
+        shape = fiona.open(inp)
+        with open(inp) as aoi:
+            aoi_resp = json.load(aoi)
+            if list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==0:
+                print('Number of current vertices '+str(len(aoi_resp['features'][0]['geometry']['coordinates'][0])))
+            elif list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==1:
+                print('Number of current vertices '+str(len(aoi_resp['features'][0]['geometry']['coordinates'][0][0])))
+            else:
+                print('Please check GeoJSON: Could not parse coordinates')
+        with fiona.Env():
+            with fiona.open(inp, 'r') as src:
+               with fiona.open(output, 'w', schema=src.schema, driver=src.driver, crs=src.crs) as sink:
+                for f in src:
+                    sink.write(vw.simplify_feature(f, number=num))
+            print('Write Completed to: '+str(output))
+    except ImportError:
+        with open(inp) as aoi:
+            aoi_resp = json.load(aoi)
+            if list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==0:
+                aoi_geom = aoi_resp['features'][0]['geometry']['coordinates']
+                print('Number of current vertices '+str(len(aoi_resp['features'][0]['geometry']['coordinates'][0])))
+            elif list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==1:
+                aoi_geom = aoi_resp['features'][0]['geometry']['coordinates'][0]
+                print('Number of current vertices '+str(len(aoi_resp['features'][0]['geometry']['coordinates'][0][0])))
+            else:
+                print('Please check GeoJSON: Could not parse coordinates')
+        ft = {
+          "type": "Feature",
+          "geometry": {
+            "type": "Polygon",
+            "coordinates": []
+          },
+          "properties": {
+            "vertex_count": num
+          }
+        }
+        ft['geometry']['coordinates']=aoi_geom
+        #returns a copy of the feature, simplified (using number of vertices)
+        b= vw.simplify_feature(ft, number=num)
+        ft['geometry']['coordinates']=b['geometry']['coordinates']
+        with open(output, 'w') as g:
+            json.dump(ft, g)
+        print('Write Completed to: '+str(output))
+#geomsimple(inp=r'C:\Users\samapriya\Downloads\vertex.geojson',output=r'C:\Users\samapriya\Downloads\v2.geojson',num=5)
diff --git a/porder/porder.py b/porder/porder.py
@@ -25,6 +25,7 @@
 import json
 import base64
 import clipboard
+from .geojson_simplify import geosimple
 from .geojson2id import idl
 from .text_split import idsplit
 from .order_now import order
@@ -58,6 +59,10 @@ def gcs_cred(cred):
 def gcs_cred_from_parser(args):
     gcs_cred(cred=args.cred)
 
+# Simplify geojson by vertex count
+def simplify_from_parser(args):
+    geosimple(inp=args.input, output=args.output, num=args.number)
+
 #Create ID List with structured JSON
 def idlist_from_parser(args):
     print('')
@@ -150,6 +155,12 @@ def main(args=None):
     required_named.add_argument('--cred', help='Path to GCS credential file', required=True)
     parser_gcs_cred.set_defaults(func=gcs_cred_from_parser)
 
+    parser_simplify = subparsers.add_parser('simplify',help='Simplifies geometry to number of vertices specified using Visvalingam-Wyatt line simplification algorithm')
+    parser_simplify.add_argument('--input',help='Input GeoJSON file')
+    parser_simplify.add_argument('--output',help='Output simplified GeoJSON file')
+    parser_simplify.add_argument('--number',help='Total number of vertices in output GeoJSON')
+    parser_simplify.set_defaults(func=simplify_from_parser)
+
     parser_idlist = subparsers.add_parser('idlist', help='Get idlist using geometry & filters')
     required_named = parser_idlist.add_argument_group('Required named arguments.')
     required_named.add_argument('--input', help='Input geometry file for now geojson/json/kml', required=True)

diff --git a/setup.py b/setup.py
@@ -20,12 +20,13 @@ def readme():
         return f.read()
 setuptools.setup(
     name='porder',
-    version='0.2.6',
+    version='0.2.7',
     packages=['porder'],
     url='https://github.com/samapriya/porder',
     package_data={'': ['bundles.json']},
     install_requires=['requests>=2.19.1','planet>=1.2.1','retrying>=1.3.3',
                       'progressbar2>=3.38.0',
+                      'visvalingamwyatt>=0.1.2',
                       'pySmartDL==1.2.5;python_version<"3.4"',
                       'pySmartDL>=1.3.1;python_version>"3.4"',
                       'shapely>=1.6.4;platform_system!="Windows"',