diff --git a/DensityProfiler_pure.py b/DensityProfiler_pure.py deleted file mode 100644 index 8da76ac..0000000 --- a/DensityProfiler_pure.py +++ /dev/null @@ -1,71 +0,0 @@ -import random -import shapefile -import pngcanvas - -def pip(x,y,poly): - n = len(poly) - inside = False - p1x,p1y = poly[0] - for i in range(n+1): - p2x,p2y = poly[i % n] - if y > min(p1y,p2y): - if y <= max(p1y,p2y): - if x <= max(p1x,p2x): - if p1y != p2y: - xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x - if p1x == p2x or x <= xinters: - inside = not inside - p1x,p1y = p2x,p2y - return inside - -# Source shapefile -r = shapefile.Reader("GIS_CensusTract_poly.shp") - -# pixel to coordinate info -xdist = r.bbox[2] - r.bbox[0] -ydist = r.bbox[3] - r.bbox[1] -iwidth = 600 -iheight = 500 -xratio = iwidth/xdist -yratio = iheight/ydist - -c = pngcanvas.PNGCanvas(iwidth,iheight,color=[255,255,255,0xff]) - -# background color -c.filledRectangle(0,0,iwidth,iheight) - -# Pen color -c.color = [139,137,137,0xff] - -# Draw the polygons -for shape in r.shapes(): - pixels = [] - for x,y in shape.points: - px = int(iwidth - ((r.bbox[2] - x) * xratio)) - py = int((r.bbox[3] - y) * yratio) - pixels.append([px,py]) - c.polyline(pixels) - -rnum = 0 -trnum = len(r.shapeRecords()) -for sr in r.shapeRecords(): - rnum += 1 - #print rnum, " of ", trnum - density = sr.record[20] - total = int(density / 50) - count = 0 - minx, miny, maxx, maxy = sr.shape.bbox - while count < total: - x = random.uniform(minx,maxx) - y = random.uniform(miny,maxy) - if pip(x,y,sr.shape.points): - count += 1 - #print " ", count, " of ", total - px = int(iwidth - ((r.bbox[2] - x) * xratio)) - py = int((r.bbox[3] - y) * yratio) - c.point(px,py,color=[255,0,0,0xff]) - -f = file("density_pure.png", "wb") -f.write(c.dump()) -f.close() - \ No newline at end of file diff --git a/PiP_Edge.py b/PiP_Edge.py deleted file mode 100644 index 360b7b9..0000000 --- a/PiP_Edge.py +++ /dev/null @@ -1,53 +0,0 @@ -# Improved point in polygon test which includes edge -# and vertex points - -def point_in_poly(x,y,poly): - - # check if point is a vertex - if (x,y) in poly: return "IN" - - # check if point is on a boundary - for i in range(len(poly)): - p1 = None - p2 = None - if i==0: - p1 = poly[0] - p2 = poly[1] - else: - p1 = poly[i-1] - p2 = poly[i] - if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]): - return "IN" - - n = len(poly) - inside = False - - p1x,p1y = poly[0] - for i in range(n+1): - p2x,p2y = poly[i % n] - if y > min(p1y,p2y): - if y <= max(p1y,p2y): - if x <= max(p1x,p2x): - if p1y != p2y: - xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x - if p1x == p2x or x <= xints: - inside = not inside - p1x,p1y = p2x,p2y - - if inside: return "IN" - else: return "OUT" - -# Test a vertex for inclusion -poligono = [(-33.416032,-70.593016), (-33.415370,-70.589604), -(-33.417340,-70.589046), (-33.417949,-70.592351), -(-33.416032,-70.593016)] -lat= -33.416032 -lon= -70.593016 - -print point_in_poly(lat, lon, poligono) - -# test a boundary point for inclusion -poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)] -x = 3 -y = 1 -print point_in_poly(x, y, poly2) diff --git a/README.md b/ProjectHome.md similarity index 84% rename from README.md rename to ProjectHome.md index 202b8da..c2217ac 100644 --- a/README.md +++ b/ProjectHome.md @@ -1,10 +1,7 @@ -# geospatialpython -Code examples from GeospatialPython.com - Sample projects and libraries from [GeospatialPython.com](http://www.GeospatialPython.com) Sample data and some scripts used in the examples on [GeospatialPython.com](http://www.GeospatialPython.com) can be found in the [Downloads](http://code.google.com/p/geospatialpython/downloads/list) section. Most sample source code is located in the [Subversion Repository](http://code.google.com/p/geospatialpython/source/browse/#svn%2Ftrunk). A good bit of the examples are based on the the [Python Shapefile Library](http://code.google.com/p/pyshp/) also hosted on Google Code. -There is a discussion group for everything above on [Google Groups](http://groups.google.com/group/geospatialpython). +There is a discussion group for everything above on [Google Groups](http://groups.google.com/group/geospatialpython). \ No newline at end of file diff --git a/PureShp2Img.py b/PureShp2Img.py deleted file mode 100644 index 7f373d2..0000000 --- a/PureShp2Img.py +++ /dev/null @@ -1,39 +0,0 @@ -""" -PureShp2Img - Convert a shapefile polygon to a PNG image using pure python - -Requires Rui Carmo's PNGCanvas module from http://the.taoofmac.com -""" - -import shapefile -import pngcanvas - -# Read in a shapefile and write png image -r = shapefile.Reader("mississippi") -xdist = r.bbox[2] - r.bbox[0] -ydist = r.bbox[3] - r.bbox[1] -iwidth = 400 -iheight = 600 -xratio = iwidth/xdist -yratio = iheight/ydist -pixels = [] -# -# Only using the first shape record -for x,y in r.shapes()[0].points: - px = int(iwidth - ((r.bbox[2] - x) * xratio)) - py = int((r.bbox[3] - y) * yratio) - pixels.append([px,py]) -c = pngcanvas.PNGCanvas(iwidth,iheight) -c.polyline(pixels) -f = file("mississippi.png","wb") -f.write(c.dump()) -f.close() -# -# Create a world file -wld = file("mississippi.pgw", "w") -wld.write("%s\n" % (xdist/iwidth)) -wld.write("0.0\n") -wld.write("0.0\n") -wld.write("-%s\n" % (ydist/iheight)) -wld.write("%s\n" % r.bbox[0]) -wld.write("%s\n" % r.bbox[3]) -wld.close \ No newline at end of file diff --git a/changeShapeType.py b/changeShapeType.py deleted file mode 100644 index 3262163..0000000 --- a/changeShapeType.py +++ /dev/null @@ -1,56 +0,0 @@ -""" -Convert one shapefile type to another -""" - -import shapefile - - -# Create a line and a multi-point -# and single point version of -# a polygon shapefile - - -## POLYLINE version -# The shapefile type we are converting to -newType = shapefile.POLYLINE - -r = shapefile.Reader("Mississippi") -w = shapefile.Writer(newType) -w._shapes.extend(r.shapes()) -# You must explicity set the shapeType of each record. -# Eventually the library will set them to the same -# as the file shape type automatically. -for s in w.shapes(): - s.shapeType = newType -w.fields = list(r.fields) -w.records.extend(r.records()) -w.save("Miss_Line") - -## MULTIPOINT version -newType = shapefile.MULTIPOINT - -w = shapefile.Writer(newType) -w._shapes.extend(r.shapes()) -for s in w.shapes(): - s.shapeType = newType -w.fields = list(r.fields) -w.records.extend(r.records()) -w.save("Miss_MPoint") - -## POINT version -newType = shapefile.POINT - -w = shapefile.Writer(newType) -# For a single point shapefile -# from another type we -# "flatten" each shape -# so each point is a new record. -# This means we must also assign -# each point a record which means -# records are usually duplicated. -for s in r.shapeRecords(): - for p in s.shape.points: - w.point(*p) - w.records.append(s.record) -w.fields = list(r.fields) -w.save("Miss_Point") \ No newline at end of file diff --git a/editshp.py b/editshp.py deleted file mode 100644 index 33c7675..0000000 --- a/editshp.py +++ /dev/null @@ -1,26 +0,0 @@ -import shapefile -# Polygon shapefile we are updating. -# We must include a file extension in -# this case because the file name -# has multiple dots and pyshp would get -# confused otherwise. -file_name = "ep202009.026_5day_pgn.shp" -# Create a shapefile reader -r = shapefile.Reader(file_name) -# Create a shapefile writer -# using the same shape type -# as our reader -w = shapefile.Writer(r.shapeType) -# Copy over the existing dbf fields -w.fields = list(r.fields) -# Copy over the existing dbf records -w.records.extend(r.records()) -# Copy over the existing polygons -w._shapes.extend(r.shapes()) -# Add a new polygon -w.poly(parts=[[[-104,24],[-104,25],[-103,25],[-103,24],[-104,24]]]) -# Add a new dbf record for our polygon making sure we include -# all of the fields in the original file (r.fields) -w.record("STANLEY","TD","091022/1500","27","21","48","ep") -# Overwrite the old shapefile or change the name and make a copy -w.save(file_name) \ No newline at end of file diff --git a/logo.png b/logo.png deleted file mode 100644 index ceed197..0000000 Binary files a/logo.png and /dev/null differ diff --git a/pointinpoly.py b/pointinpoly.py deleted file mode 100644 index f6c018a..0000000 --- a/pointinpoly.py +++ /dev/null @@ -1,34 +0,0 @@ - -# Determine if a point is inside a given polygon or not -# Polygon is a list of (x,y) pairs. This fuction -# returns True or False. The algorithm is called -# "Ray Casting Method". - -def point_in_poly(x,y,poly): - - n = len(poly) - inside = False - - p1x,p1y = poly[0] - for i in range(n+1): - p2x,p2y = poly[i % n] - if y > min(p1y,p2y): - if y <= max(p1y,p2y): - if x <= max(p1x,p2x): - if p1y != p2y: - xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x - if p1x == p2x or x <= xinters: - inside = not inside - p1x,p1y = p2x,p2y - - return inside - -## Test - -polygon = [(0,10),(10,10),(10,0),(0,0)] - -point_x = 5 -point_y = 5 - -## Call the fuction with the points and the polygon -print point_in_poly(point_x,point_y,polygon) \ No newline at end of file diff --git a/pyshpLargeWriter.py b/pyshpLargeWriter.py deleted file mode 100644 index ac30adf..0000000 --- a/pyshpLargeWriter.py +++ /dev/null @@ -1,186 +0,0 @@ -import cStringIO -import struct -import shapefile -import random -import time - -class LargeWriter: - def __init__(self, filename=None, shapeType=1, hasShx=True): - self.filename = filename - self.hasShx = hasShx - # Count records for metadata - self.count = 0 - # Maximum number of records before disk flush - self.max = 1000 - self.started = False - self.minx = 0 - self.miny = 0 - self.maxx = 0 - self.maxy = 0 - self.numRecs = 0 - self.tmpShp = cStringIO.StringIO() - if self.hasShx: - self.tmpShx = cStringIO.StringIO() - self.tmpDbf = cStringIO.StringIO() - self.shp = open("%s.shp" % self.filename, "wb") - if self.hasShx: - self.shx = open("%s.shx" % self.filename, "wb") - self.dbf = open("%s.dbf" % self.filename, "wb") - self.dbfHdrLen = 0 - self.w = shapefile.Writer(shapeType) - - def endcap(self): - """First and last batches""" - self.started = True - self.tmpShp.seek(36) - self.xmin = struct.unpack("i", self.tmpShx.read(4))[0] - newOffset = adder + offset - self.tmpShx.seek(-4, 1) - self.tmpShx.write(struct.pack(">i", newOffset)) - self.tmpShx.seek(8, 1) - shxcur = self.tmpShx.tell() - # Get shp bounding box - self.tmpShp.seek(36) - xmin = struct.unpack("= self.max: - if self.hasShx: - self.w.save(shp=self.tmpShp, shx=self.tmpShx, dbf=self.tmpDbf) - else: - self.w.save(shp=self.tmpShp, dbf=self.tmpDbf) - if not self.started: - self.endcap() - else: - self.batch() - # Reset the buffers for the next batch - self.tmpShp = cStringIO.StringIO() - if self.hasShx: - self.tmpShx = cStringIO.StringIO() - self.tmpDbf = cStringIO.StringIO() - self.count = 0 - - def save(self): - self.shp.seek(0,2) - shpLength = self.shp.tell() / 2 - self.shp.seek(24) - self.shp.write(struct.pack(">i", shpLength)) - self.shp.seek(36) - self.shp.write(struct.pack("<4d", self.xmin,self.ymin,self.xmax,self.ymax)) - if self.hasShx: - self.shx.seek(0,2) - shxLength = self.shx.tell() / 2 - self.shx.seek(24) - self.shx.write(struct.pack(">i", shxLength)) - self.shx.seek(36) - self.shx.write(struct.pack("<4d", self.xmin,self.ymin,self.xmax,self.ymax)) - # update dbf record count - self.dbf.seek(4) - self.dbf.write(struct.pack(" last: - print "%s %% done - Shape %s/%s at %s" % (percent, count, total, time.asctime()) - return percent - - -# Create a large shapefile writer with -# a filname and a shapefile type. -# 1=Point, 5=Polygon -lw = LargeWriter("giantShp", 1, hasShx=True) -#lw = LargeWriter("giantShp", 1, hasShx=False) - -# Add some dbf fields -lw.w.field("ID", "C", "40") -lw.w.field("X", "C", "40") -lw.w.field("Y", "C", "40") - -# Progress counter -status = 0 - -# Number of random points to write -total = 107374177 - -for i in range(total): - # Progress meter - status = progress(status, i, total) - # Generate a random point - x,y = random_point() - x1,y1 = random_point() - x2,y2 = random_point() - x3,y3 = random_point() - x4,y4 = random_point() - # Call the shapefile point() method - lw.w.point(x,y) - # Call the specialized record method - lw.record(i,x,y) - #lw.w.poly(parts=[[[x,y],[x1,y1],[x2,y2],[x3,y3],[x4,y4]]]) - #lw.record(i,x,y) - #x,y = random_point() - #x1,y1 = random_point() - #x2,y2 = random_point() - #x3,y3 = random_point() - #x4,y4 = random_point() - #Call the shapefile poly() method - #lw.w.poly(parts=[[[x,y],[x1,y1],[x2,y2],[x3,y3],[x4,y4]]]) - #Call the specialized record method - #lw.record(i,x,y) - -# Special LargeWriter save method -lw.save() \ No newline at end of file diff --git a/shp2img.py b/shp2img.py deleted file mode 100644 index a355f60..0000000 --- a/shp2img.py +++ /dev/null @@ -1,41 +0,0 @@ -""" -shp2img.py - creates a png image and world file (.pgw) from a shapefile -containing a single polygon. -author: jlawhead@geospatialpyton.com -date: 20101204 - -This script requires the Python Imaging Library. -The sample shapefile used is available at: -http://geospatialpython.googlecode.com/files/Mississippi.zip -""" - -import shapefile -import Image, ImageDraw - -# Read in a shapefile -r = shapefile.Reader("mississippi") -xdist = r.bbox[2] - r.bbox[0] -ydist = r.bbox[3] - r.bbox[1] -iwidth = 400 -iheight = 600 -xratio = iwidth/xdist -yratio = iheight/ydist -pixels = [] -for x,y in r.shapes()[0].points: - px = int(iwidth - ((r.bbox[2] - x) * xratio)) - py = int((r.bbox[3] - y) * yratio) - pixels.append((px,py)) -img = Image.new("RGB", (iwidth, iheight), "white") -draw = ImageDraw.Draw(img) -draw.polygon(pixels, outline="rgb(203, 196, 190)", fill="rgb(198, 204, 189)") -img.save("mississippi.png") - -# Create a world file -wld = file("mississippi.pgw", "w") -wld.write("%s\n" % (xdist/iwidth)) -wld.write("0.0\n") -wld.write("0.0\n") -wld.write("-%s\n" % (ydist/iheight)) -wld.write("%s\n" % r.bbox[0]) -wld.write("%s\n" % r.bbox[3]) -wld.close \ No newline at end of file diff --git a/subset.py b/subset.py deleted file mode 100644 index 873aacd..0000000 --- a/subset.py +++ /dev/null @@ -1,20 +0,0 @@ -import shapefile - -# Create a reader instance -r = shapefile.Reader("buildingfootprints/Building_Footprint") -# Create a writer instance -w = shapefile.Writer(shapeType=shapefile.POLYGON) -# Copy the fields to the writer -w.fields = list(r.fields) -# Grab the geometry and records from all features -# with the correct county name -selection = [] -for rec in enumerate(r.records()): - if rec[1][1].startswith("Hancock"): - selection.append(rec) -# Add the geometry and records to the writer -for rec in selection: - w._shapes.append(r.shape(rec[0])) - w.records.append(rec[1]) -# Save the new shapefile -w.save("buildingfootprints/HancockFootprints") diff --git a/zipshape.py b/zipshape.py deleted file mode 100644 index ab41007..0000000 --- a/zipshape.py +++ /dev/null @@ -1,45 +0,0 @@ -""" -Example of saving a shapefile to a file-like -object using a zip file as the target. Uses -the individual save methods for each type -of file instead of the general pyshp -save method. -""" - -import zipfile -import StringIO -import shapefile - -# Set up buffers for saving -shp = StringIO.StringIO() -shx = StringIO.StringIO() -dbf = StringIO.StringIO() - -# Make a point shapefile -w = shapefile.Writer(shapefile.POINT) -w.point(90.3, 30) -w.point(92, 40) -w.point(-122.4, 30) -w.point(-90, 35.1) -w.field('FIRST_FLD') -w.field('SECOND_FLD','C','40') -w.record('First','Point') -w.record('Second','Point') -w.record('Third','Point') -w.record('Fourth','Point') - -# Save shapefile components to buffers -w.saveShp(shp) -w.saveShx(shx) -w.saveDbf(dbf) - -# Save shapefile buffers to zip file -# Note: zlib must be available for -# ZIP_DEFLATED to compress. Otherwise -# just use ZIP_STORED. -z = zipfile.ZipFile("myshape.zip", "w", zipfile.ZIP_DEFLATED) -z.writestr("myshape.shp", shp.getvalue()) -z.writestr("myshape.shx", shx.getvalue()) -z.writestr("myshape.dbf", dbf.getvalue()) -z.close() -