HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

Python program to check if a set of points is at land or at sea

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
seapointsprogramlandpythoncheckset

Problem

I have written a program that divides an geographical area up in cells with size 0.002 degrees longitude by 0.001 degrees latitude, and check every cell whether it's a land area or a sea area. I've done this to effectively to get all the sea area of that area. This is probably insane from a GIS standpoint. Please review my code, call me an insane idiot and help me to get it better.

"""Function to iterate with decimal step size"""
def seq(start, end, step):
    assert(step != 0)
    sample_count = abs(end - start) / step
    return itertools.islice(itertools.count(start, step), sample_count)

if __name__ == "__main__":
    bottomlat = 70.390829
    toplat = 70.855549
    bottomlong = 27.692793
    toplong =  28.728254
    gridsize = 0.001
    map = Basemap(
    projection="merc",
    resolution = 'h', area_thresh = 0.001,
    llcrnrlon=bottomlong, llcrnrlat=bottomlat,
    urcrnrlon=toplong, urcrnrlat=toplat)
    map.drawcoastlines(color='black')
    """Divide global area into grid cells"""
    for i in seq(bottomlat, toplat, gridsize):
        for j in seq(bottomlong, toplong, gridsize*2):
            lon, lat = map(j,i)
            if not map.is_land(lon,lat):
                plt.scatter(lon,lat, s=0.1, color='c', edgecolors='none')
    plt.show()


The ultimate use case is:
I've god millions of lat/long positions from ships. I want to see how large percentage of the sea that has a lat/long position from a ship (cells with a ship position divided by total cell of ocean). So my initial thought was to check each cell if its sea, then perform a check if there is a ship position there.

Solution

Style

Your docstring at the top of the seq function should be moved underneath the function signature, like this:

def seq(start, end, step):
    """Function to iterate with decimal step size"""
    ...


The variable map, should be written like this, to improve clarity and readability:

map = Basemap(
    projection="merc",
    resolution = 'h', 
    area_thresh = 0.001,
    llcrnrlon=bottomlong, 
    llcrnrlat=bottomlat,
    urcrnrlon=toplong, 
    urcrnrlat=toplat
)


All these variables:

bottomlat = 70.390829
toplat = 70.855549
bottomlong = 27.692793
toplong =  28.728254
gridsize = 0.001


Should be renamed so they appropriately display the fact that they are constants.

BOTTOM_LAT = 70.390829
TOP_LAT = 70.855549
BOTTOM_LONG = 27.692793
TOP_LONG =  28.728254
GRID_SIZE = 0.001


This comment:

"""Divide global area into grid cells"""


Should be changed into an inline comment:

# Divide global area into grid cells


For more help with style, visit PEP8, Python's official style guide.

Error handling in seq

You shouldn't be using assert to check if a variable meets certain constraints unless you're doing testing. The proper way would be to do something like this:

if step != 0:
    sample_count = abs(end - start) / step
    return itertools.islice(itertools.count(start, step), sample_count)
raise ValueError("Step must not be equal to zero.")

Code Snippets

def seq(start, end, step):
    """Function to iterate with decimal step size"""
    ...
map = Basemap(
    projection="merc",
    resolution = 'h', 
    area_thresh = 0.001,
    llcrnrlon=bottomlong, 
    llcrnrlat=bottomlat,
    urcrnrlon=toplong, 
    urcrnrlat=toplat
)
bottomlat = 70.390829
toplat = 70.855549
bottomlong = 27.692793
toplong =  28.728254
gridsize = 0.001
BOTTOM_LAT = 70.390829
TOP_LAT = 70.855549
BOTTOM_LONG = 27.692793
TOP_LONG =  28.728254
GRID_SIZE = 0.001
"""Divide global area into grid cells"""

Context

StackExchange Code Review Q#97748, answer score: 6

Revisions (0)

No revisions yet.