GeospatialPython.com

Can you advise as to how to convert the GPS time o...

2012-05-07T10:19:04.712-05:00

Can you advise as to how to convert the GPS time of week to a complete time structure like SYSTEMTIME? Seems like the GPS week number would be needed for this but it does not seem to be recorded. Thanks!

There is a standard method that works every time. ...

2012-05-03T12:10:42.426-05:00

There is a standard method that works every time.

You the visual studio solution downloadable from http://www.inferentialpower.com/Bio/Ing/ObjectParser/.

It contains an "InPolygon" method written in C#.

The solution also contains an "Imprint" method for imprinting polygons onto bitmaps.

bluejp - Sorry for the slow reply. That is really ...

2012-04-13T10:53:36.504-05:00

bluejp - Sorry for the slow reply. That is really weird. What platform are you on? When you call Editor.save() you are calling the Writer.save() method which explicitly closes all file handles.

All the Editor class does is wrap the Reader and Writer classes so you don't have to explicitly create separate reader and writer objects to change a shapefile. But of course you can do that manually. You just create a reader and a writer and copy everything into the writer and make changes, then save over the old file. There are several posts where I do this.

Also - the Editor is not as well tested so people continually find bugs but it's slowly getting better. The Reader and Writer classes are pretty solid and should work. I have examples of this technique in the following posts:
http://geospatialpython.com/2011/02/merging-lots-of-shapefiles-quickly.html

http://geospatialpython.com/2010/12/subsetting-shapefile-by-attributes.html

I would try that and see if you get the same results. If you're using it in a web application it could be some web server behavior. But I really don't know. Let me know what you figure out and I'll see if I can help troubleshoot.

This is great write-up. I want to find it out so I...

2012-04-11T11:53:37.103-05:00

This is great write-up. I want to find it out so I'm looking to download it to the link you gave. Thanks.

Accounting Packages

Joel... great job you have done here. Question.. w...

2012-04-10T11:01:05.804-05:00

Joel... great job you have done here. Question.. when working with the Editor... is there any function to close the editing operations. I experience that the Editor "locks" the files, even after closing python. Any suggestions?

I would like to use the code to clip a 1-band rast...

2012-03-22T05:05:56.920-05:00

I would like to use the code to clip a 1-band raster file based on a multi-features polygon shapefile, but some errors are happening, like:

clip = srcArray[:, ulY:lrY, ulX:lrX]
IndexError: too many indices

VeeAnn - I'm looking into this issue. It may ...

2012-03-14T14:56:25.122-05:00

VeeAnn - I'm looking into this issue. It may be a bug.

I would also like to know if this is possible.

2012-03-14T13:37:32.536-05:00

I would also like to know if this is possible.

I created what appears to be a polylineM shapefile...

2012-03-09T10:13:41.721-06:00

I created what appears to be a polylineM shapefile. It doesn't seem quite right in that none of the measurements seem to have gone over, but nonetheless, I can open with Esri and it thinks it's a polylineM. However, I can't seem to read the records using this method. So my problem is two-fold. I need a legitimate polylineM and I need to be able to read it. I'm sure I'm doing something wrong, but would appreciate any assistance.

I think I have the second algorithm figured out. I...

2012-03-03T09:46:42.995-06:00

I think I have the second algorithm figured out. I can successfully bin my 3000 test shapefiles (randomly generated with fixed extents). I'm still getting some random misbins on the World Cities shapefile but I think these are due to minor errors in implementation.

In a nutshell, the algorithm examines a bin and checks to see if its grandchildren contain features. If not, it pulls up any features from its child nodes.

Francisco, Yes, that is basically the definition ...

2012-03-02T09:14:58.934-06:00

Francisco,

Yes, that is basically the definition of binary division and we have been using it for a while. You can see it in this photo which shows the 1 integer unit seams between cells.

However, if you test enough shapefiles with it you will start to find a small percentage of seam features do not end up where you expect them to be. Typically, they will be one to two levels higher in the tree. There appears to be a second (in the sense of additional, not different) algorithm that is run on the seam features, probably for some sort of optimization.

Mark, The seam problem is not a second algorithm....

2012-03-01T21:31:34.890-06:00

Mark,

The seam problem is not a second algorithm. I figure it out back when I was trying to understand what's going on. It is not exactly a continous split, but rather a discrete one. Following these instructions you should be able to recreate the tree with no seam problem.

Suppose we are doing the first vertical split, thus we have the whole shapefile bounding box, with xmin=0 and xmax=255. Then you calculate center = xmin + (xmax - xmin) / 2 + 1, where the division is an integer division, not real. As you shall see in a bit, this will always result in an even number.

Suppose we have a shape and that shape.xmin and shape.xmax are the byte range for the x's. Then we have shape.xmin < shape.xmax, even if the feature is just a point. This is because when changing the real number to a byte, the xmin gets rounded down and xmax gets rounded up.

To place the shape in the proper node in the tree, do the following. If shape.xmin <= center and shape.xmax>center, then the shape bounding box is assumed to touch the vertical line, even if the not-rounded bounding box does not touch the exactly vertical line. If shape.xmin > center, then the shape is placed on the first child or one of its descendant, and for the next vertical split xmin=center. If shape.xmax <= center, then the shape is placed on the second child and for the next vertical split xmax=center-1.

Suppose we start with xmin=0 and xmax=255. Then center=128. Say we now want to search on the right. Then xmin=128 and xmax=255 and center=192. Say we want to go to the left now. Then xmin=128 and xmax=191 and center=160. xmin will always be an even number and xmax an odd number, in fact the width (xmax - xmin) is always a power of 2 minus one (255, 127, 63, 31, etc), which is what causes the center to always be even.

The same principle applies to horizontal splits. I followed this "discrete" splitting scheme and checked with all the features in my shapefiles, and I got them all placed in the correct bins. Please let me know if it works for you too.

Awesome! I was hoping you'd add something like...

2012-03-01T11:10:01.498-06:00

Awesome! I was hoping you'd add something like this. I just started using Pyshp in the past 3 weeks.

Francisco, Yes, it is a binary division. It's...

2012-03-01T11:09:29.539-06:00

Francisco,

Yes, it is a binary division. It's essentially an "Interval KD Tree" but the splits are at the midpoint of the node's range, not at the median point in the set. Also, there's no tree applied to the seam features.

There is a splitting threshold of 8 features per cell initially. However, the tree size is fixed based on the total number of features. It selects the size of the tree so that there are 8 features or less per cell on average.

Because you can only split the range of the 8 bit mapping so many times, the maximum size of the tree is fixed. How the space is split is also a bit odd. It doesn't split on a single value - there's actually a 1 unit "gutter" between the cells. There's also padding on the cells on the low edges.

Once you implement all this you'll run into some strange behavior - features on the seams between splits won't end up in the bins where your algorithm predicts. It appears there's a second algorithm for the seam features. I've figured it out, but it only seems to apply to first 6 levels of the tree. I can see the justification for why there could be two algorithms, but I haven't completely solved the problem.

But we are very close.

2012-03-01T10:04:40.830-06:00

This comment has been removed by the author.

First of all thanks for what you have done so far ...

2012-03-01T10:04:14.369-06:00

First of all thanks for what you have done so far in the reverse ingineering process of sbn file. It help me a lot with the C# library I built to read shapefiles. Here are my two cents to this process

The sbn is an spatial index based on a binary tree, where each node contains one or more boxes. When a node has too many boxes it splits. The split draws a line in the middle. Any box that touches the line is kept in the node, while any box on one side of the line is sent to the first child, and any box on the other side of the line is sent to the second child. A vertical split divides the square with a vertical line into two rectangles, and its first child receives the boxes to the right of the line and the second child receives the boxes to the left of the line. A horizontal split divides the rectangle with a horizontal line into two squares. The boxes above the line are sent to the first child while the boxes below to the second child. The root is a vertical split. The children of a vertical split are always horizontal splits, while the children of a horizontal split are always vertical splits.

The list of bins is an array contaning the tree. the first element is the root node, the second and third elements its children the 4th and 5th are the children of the the 2nd bin and the 6th and 7th element are the children of the third bin, and so on. If a node has not been split, then all its descendants will be empty bin records. If a node has been split, but no box touches the splitting line, then than bin will be an empty bin record.

The one thing I have not figure out yet is what's the splitting threshold, i.e., when a node is considered full and splitting occur. The files I have looked at have different numbers. One of them, which had a few shapes, does not have more than 30 boxes in a node, while another that has lots of shapes has no more than 700 shapes in a node. This is a problem when writing the sbn file. However is not a problem when reading it and using it. With the large file, when I read it without the index it took a few seconds to retrive the shapes that intersected certain small box, but when I started using the index file, it worked in the blink of an eye.

Very nice piece of code and very useful to merge a...

2012-02-22T03:06:33.966-06:00

Very nice piece of code and very useful to merge a lot of tiled data.

I however had to update a little bit the shapefile.py (date: 20110927, version: 1.1.4) to process my pointZ shapefiles. On line 699, I've change s.points[0][2] by s.z[0] and on line 705 s.points[0][3] by s.m. Does it seems correct?

Thanks a lot,
David

Mike, There's no conversion or measurement go...

2012-02-08T09:15:01.496-06:00

Mike,

There's no conversion or measurement going on so I think this is fine for lat/long. It should be abstract enough and works on arbitrary polygons so I think it would be ok. Correct me if you discover otherwise.

Joel, I would be happy to provide my script, but ...

2012-02-07T22:27:39.068-06:00

Joel,

I would be happy to provide my script, but I want to make sure I have my head on straight. Don't you need to convert the lat/lon points into points on a Cartesian coordinate system for this point-in-polygon algorithm to work?

Antonio, You shouldn't have to extract the po...

2012-02-06T14:30:09.218-06:00

Antonio,

You shouldn't have to extract the polygons to separate files.

This script could be modified to do that. Instead of "poly = lyr.GetNextFeature()", you would build and array of polygons for each feature in the shapefile.

And then you would need to put the "Map Points to Pixels" section in a loop to apply each polygon to the clipping mask.

I haven't tried this trick but it should work.

When I have more time I can try to create an example. If you get it working let me know.

Mike, Sorry for the slow reply. I don't have...

2012-02-06T14:13:00.270-06:00

Mike,

Sorry for the slow reply. I don't have anything coded up but here's some pretty straight forward math that may help:

http://tinyurl.com/7cqzwdw

If you implement this in Python please let me know. I'd love to share that with the community.

Joel: Do you have any suggestions for determining ...

2012-02-03T20:00:53.875-06:00

Joel: Do you have any suggestions for determining (with Python) if a convex or concave polygon runs clockwise? Your script works well, but the set of polygons I have are a mixture of clockwise and counterclockwise. Thank you in advance.

Thank you for the post, If I have a single shapef...

2012-02-03T12:04:23.428-06:00

Thank you for the post,

If I have a single shapefile with several polygons inside it how can I make several clips for each of that polyon? I have to extract them from a shapefile and iterate over them?

Does anyone know how to do this in R?

2012-01-25T11:10:12.825-06:00

Does anyone know how to do this in R?

oops, not the shp file, the error said the dbf is ...

2012-01-17T10:39:45.826-06:00

oops, not the shp file, the error said the dbf is corrupt or does not have same number of objects as the shx file, run your program and it works, pretty simple, I wonder why its not working properly though, any ideas?