Python program that helps its user explore a polygon shapefile. The program should first ask the user to provide the path to a polygon shapefile. If the shapefile does not contain polygons or multipolygon, the program should tell the user to provide a valid polygon shapefile and then exit. When the file is valid, the program provides some basic information about the file, including

• The number of features,
• The number of attributes (properties), and
• The name of each attribute.

Then the program goes on a loop to ask the user to input an index number of a polygon feature. If a valid number is entered, the program should print the following info:

• The name and value for all attributes,
• Type of polygon (Polygon or MultiPolygon), and
• The total number of points of all parts and rings.

If an invalid index is entered, the program should print "Invalid index" and then ask the index again (in the loop). The loop ends when the user enters the word DONE or blank (empty string). Before the program ends, it should print out something like "Thank you for using our program!"

To get you started, the following is the beginning of the program that asks for the valid polygon shapefile. You can comment out the line of using input and uncomment the line underneath it (with the correct path to the right shapefile on your own computer). In this way, you don't have to enter the shapefile every time you run the code, which will save you tremendous time when testing your code. Just reset those comments when you are done before you submit the final code.

import sys sys.path.append('/home/xiao/lib/gisalgs') from geom.shapex import * is_polygon_shp = False while not is_polygon_shp: fname = input("Enter a shapefile name: ") # fname = '/home/xiao/lib/gisalgs/data/uscnty48area.shp' try: mapdata = shapex(fname) except Exception as e: print('Error:', e) print('Make sure your enter a valid shapefile') continue shp_type = mapdata[0]['geometry']['type'] if shp_type in ['Polygon', 'MultiPolygon']: is_polygon_shp = True else: print('Make sure it is a polygon or multipolygon shapefile') print()

In your answer, please present both the code and plenty of examples of how your program runs on invalid and valid inputs.

There is no need to include the data (the shapefile you used). To provide examples of how your program works, you can simply copy and past the process of using the program as Python comments (using triple quotes). Below is a sample session of running of the program.

Enter a shapefile name: wrong name! Error: Need a .shp file. Make sure your enter a valid shapefile Enter a shapefile name: /home/xiao/lib/gisalgs/data/ne_110m_coastline.shp Make sure it is a polygon or multipolygon shapefile Enter a shapefile name: /home/xiao/lib/gisalgs/data/uscnty48area.shp There are 3109 features in the shapefile. There are 8 fieds. The ields are: NAME STATE_NAME FIPS UrbanPop Area AreaKM2 GEO_ID PopDensity Enter feature index: 0 Properties: {'NAME': 'Gallatin', 'STATE_NAME': 'Montana', 'FIPS': '30031', 'UrbanPop': 56.73512111, 'Area': 6817105293.22, 'AreaKM2': 6817.10529322, 'GEO_ID': '05000US30031', 'PopDensity': 9}. This is a Polygon Number of rings is 1 and total number of points is 72 Enter feature index: -1 Invalid index: need an integer Enter feature index: 1000 Properties: {'NAME': 'Kearny', 'STATE_NAME': 'Kansas', 'FIPS': '20093', 'UrbanPop': 0.0, 'Area': 2255830959.72, 'AreaKM2': 2255.83095972, 'GEO_ID': '05000US20093', 'PopDensity': 2}. This is a Polygon Number of rings is 1 and total number of points is 10 Enter feature index: 2788 Properties: {'NAME': 'Rockbridge', 'STATE_NAME': 'Virginia', 'FIPS': '51163', 'UrbanPop': 3.551518647, 'Area': 1556498272.31, 'AreaKM2': 1556.49827231, 'GEO_ID': '05000US51163', 'PopDensity': 13}. This is a Polygon Number of rings is 3 and total number of points is 51 Enter feature index: 2671 Properties: {'NAME': 'Suffolk', 'STATE_NAME': 'New York', 'FIPS': '36103', 'UrbanPop': 96.98274374, 'Area': 6146392494.3, 'AreaKM2': 6146.3924943, 'GEO_ID': '05000US36103', 'PopDensity': 230}. This is a MultiPolygon Number of parts is 2 and total number of points is 35 Enter feature index: huh? Invalid index: need an integer Enter feature index: 100000 Invalid index: out of range Enter feature index: done Invalid index: need an integer Enter feature index: done? Invalid index: need an integer Enter feature index: 90 Properties: {'NAME': 'Toole', 'STATE_NAME': 'Montana', 'FIPS': '30101', 'UrbanPop': 57.43307386, 'Area': 5039524604.96, 'AreaKM2': 5039.52460496, 'GEO_ID': '05000US30101', 'PopDensity': 1}. This is a Polygon Number of rings is 1 and total number of points is 21 Enter feature index: DONE Thank you for using our program!

shapex.py:

''' A class and related functions that handle reading of shapefiles It now only supports four types of shapefile: Point, MultiPoint (not tested), PolyLine, Polygon. History October 9, 2017 Support slicing! October 8, 2017 First version. Supports four types. No slicing Credit: the read_dbf method is adopted from http://code.activestate.com/recipes/362715/ Author Ningchuan Xiao n..o@gmail.com ''' from struct import unpack, calcsize from os.path import isfile from datetime import date shapefile_types = { 0: 'Null Shape', 1: 'Point', 3: 'PolyLine', 5: 'Polygon', 8: 'MultiPoint', 11: 'PointZ', 13: 'PolyLineZ', 15: 'PolygonZ', 18: 'MultiPointZ', 21: 'PointM', 23: 'PolyLineM', 25: 'PolygonM', 28: 'MultiPointM', 31: 'MultiPatch' } supported_types = [ 'Point', 'MultiPoint', 'PolyLine', 'Polygon' ] def clockwise(polygon): '''calculate 2*A polygon: [ [x, y], [x, y], ... ] polygon = [ [1, 0], [2,0], [2,2], [1,2], [1, 0] ] clockwise(polygon) # False polygon.reverse() clockwise(polygon) # True ''' if polygon[0] != polygon[-1]: return num_point = len(polygon) A = 0 for i in range(num_point-1): p1 = polygon[i] p2 = polygon[i+1] ai = p1[0] * p2[1] - p2[0] * p1[1] A += ai return A<0 class shapex: ''' A class for points in Cartesian coordinate systems. Examples >>> fname = '/Users/xiao/lib/gisalgs/data/uscnty48area.shp' >>> shp = shapex(fname) >>> print(shp[60]) ''' def __init__(self, fname): if not fname.endswith('.shp'): raise Exception('Need a .shp file.') self.fname_shp = fname self.fname_shx = fname[:-3]+'shx' self.fname_dbf = fname[:-3]+'dbf' if not isfile(self.fname_shp) or not isfile(self.fname_shx) or not isfile(self.fname_dbf): raise Exception('Need at least three files: .shp, .shx, .dbf') self.open_shapefile() def open_shapefile(self): self.f_shx = open(self.fname_shx, 'rb') h1 = unpack('>7i', self.f_shx.read(28)) h2 = unpack('<2i 8d', self.f_shx.read(72)) file_length = h1[-1] self.num_rec = (file_length-50)//4 self.f_shp = open(self.fname_shp, 'rb') h1 = unpack('>7i', self.f_shp.read(28)) # BIG h2 = unpack('<2i 8d', self.f_shp.read(72)) # LITTLE self.file_length = h1[-1] self.version = h2[0] self.shape_type = shapefile_types[h2[1]] # self.xmin, self.ymin, self.xmax, self.ymax, self.zmin, self.zmax, self.mmin, self.mmax = h2[2:10] self.xmin = h2[2] self.ymin = h2[3] self.xmax = h2[4] self.ymax = h2[5] self.zmin = h2[6] self.zmax = h2[7] self.mmin = h2[8] self.mmax = h2[9] self.this_feature_num = 0 # get (offset, content length) pairs from shx # remember each record has a header of 8 bytes index = unpack('>'+'i'*self.num_rec*2, self.f_shx.read(self.num_rec*4*4)) self.index = [(index[i]*2, index[i+1]*2) for i in range(0, len(index), 2)] # get schema, etc. self.f_dbf = open(self.fname_dbf, 'rb') dbf_numrec, lenheader = unpack('self.num_rec: raise Exception('Feature index out of range (' + str(i) + ')') pos = self.index[i] self.f_shp.seek(pos[0] + 8) # skip record hearder, which is not useful if self.shape_type == 'Polygon': feature = self.readpolygon() if self.shape_type == 'PolyLine': feature = self.readpolygon() if feature['geometry']['type'] == 'MultiPolygon': feature['geometry']['type'] = 'MultiLineString' else: feature['geometry']['type'] = 'LineString' if self.shape_type == 'Point': feature = self.readpoint() if self.shape_type == 'MultiPoint': feature = self.readmultipoint() # get properties here. properties = self.read_dbf(i) feature['properties'] = properties feature['id'] = i return feature else: raise TypeError('Invalid index') def robust_decode(self, bs): ''' https://stackoverflow.com/questions/24475393/unicodedecodeerror-ascii-codec-cant-decode-byte-0xc3-in-position-23-ordinal Convert a byte string to unicode. Try UTF8 first, if not working then latin1. ''' cr = None try: cr = bs.decode('utf8') except UnicodeDecodeError: cr = bs.decode('latin1') return cr def read_dbf(self, i): # Note: dtypes of D, L are note tested self.f_dbf.seek(self.dbf_header_length + i * self.formatsize) record = unpack(self.formatstr, self.f_dbf.read(self.formatsize)) if record[0] == ' ': return ' ' * self.formatsize result = [] for (name, dtype, size, deci), value in zip(self.fields, record): value = value.decode('latin1') # value = value.decode('ascii') // works for Python 2 # value = self.robust_decode(value) if name == 'DeletionFlag': continue if dtype == 'N': value = value.replace('\0', '').lstrip() if value == '': value = 0 elif deci: value = float(value) else: value = int(value) elif dtype == 'C': value = value.rstrip() elif dtype == 'D': y, m, d = int(value[:4]), int(value[4:6]), int(value[6:8]) value = date(y, m, d) elif dtype == 'L': value = (value in 'YyTt' and 'T') or (value in 'NnFf' and 'F') or '?' elif dtype == 'F': value = float(value) result.append(value) properties = {} for fi in range(1, len(self.fields)): properties[self.fields[fi][0]] = result[fi-1] return properties def readpoint(self): point = unpack('= self.num_rec: self.this_feature_num = 0 raise StopIteration feature = self.__getitem__(self.this_feature_num) self.this_feature_num += 1 return feature def close(self): self.f_shp.close() self.f_shx.close() self.f_dbf.close() @property def bounds(self): return(self.xmin, self.ymin, self.xmax, self.ymax) @property def schema(self): myschema = {} myschema['geometry'] = self.shape_type properties = [] for fi in range(1, len(self.fields)): name = self.fields[fi][0] f1 = self.fields[fi][1] f2 = self.fields[fi][2] dci = self.fields[fi][3] if f1 == 'C': fmt = 'str:' + str(f2) elif f1 == 'F': fmt = 'float:' + str(f2) + '.' + str(dci) elif f1 == 'N': if dci == 0: fmt = 'int:' + str(f2) else: fmt = 'float:' + str(f2) + '.' + str(dci) elif f1 == 'D': fmt = 'datetime' else: fmt = 'other' properties.append((name, fmt)) myschema['properties'] = properties return myschema if __name__ == '__main__': fname = '/Users/printlnchui/Desktop/GEOG_5222/uscnty48area.shp' # fname = '/Users/xiao/lib/gisalgs/data/ne_110m_populated_places.shp' # fname = '/Users/xiao/lib/gisalgs/data/ne_110m_admin_0_countries.shp' shp = shapex(fname) print('Number of fetures:', len(shp)) # this tests all the geometry types types = {} for f in shp: t = f['geometry']['type'] if t[:5] != 'Multi': if len(f['geometry']['coordinates']) > 1: t = t + '_Parts' if t not in types: types[t] = 1 else: types[t] += 1 print(types) print('Shape type:', shp.shape_type) print('Schema: ', shp.schema) print('Bunds: ', shp.bounds) shp.close()

Here is the link for the project: https://www.dropbox.com/s/d0mlcva32bwye0q/shapex_geojson.html?dl=0