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

Question

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.

Answer 1

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.")

---

Answer 2

Superficial issues

"""Function to iterate…""" is a docstring, which belongs inside the function definition, not before it.

"""Divide global area into grid cells""" should be written as a comment (# Divide global area into grid cells).

The way you formatted map = Basemap(…) makes it hard to see that it's one statement over five lines. The four continuation lines should be indented some more. PEP 8 says that resolution='h' and area_thresh=0.001 should not have spaces around the equals signs.

Naming

map unfortunately shadows the built-in function with the same name, but I can live with that.

"Longitude" is sometimes shortened to …lon, and sometimes to …long.

I find bottomlong and toplong confusing, because they are actually left-right boundaries.

lon, lat = map(j, i) is very misleading: those are points in the map coordinate system, not longitude and latitude! I would consider i and j to be poor names as well, since they have the connotation of being integers. The loop should therefore be written as:

# Divide global area into grid cells
for lon in seq(minlon, maxlon, 2 * gridsize):
    for lat in seq(minlat, maxlat, gridsize):
        x, y = map(lon, lat)
        if not map.is_land(x, y):
            plt.scatter(x, y, s=0.1, color='c', edgecolors='none')

Just the map

If your goal is just to obtain a nicely colored map of Tanafjord, just use the color-filling functions.

map.drawcoastlines()
map.drawmapboundary(fill_color='c')
map.fillcontinents(color='white', lake_color='c')

Sampling

On the other hand, if you are interested in finding the land-vs.-sea ratio using sampling, the most relevant function would be maskoceans(). Instead of defining seq(), you should use makegrid(). You do not need to write any loops at all.

from mpl_toolkits.basemap import Basemap, maskoceans
import matplotlib.pyplot as plt
import numpy as np

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
)

lons, lats, x, y = map.makegrid(int((toplong - bottomlong) / gridsize / 2),
                                int((toplat - bottomlat) / gridsize),
                                returnxy=True)

# Estimate area of sea by sampling
samples = len(lons) * len(lons[0])
ocean = maskoceans(lons, lats, datain=np.arange(samples),
                   resolution='h', grid=1.25)
ocean_samples = np.ma.count_masked(ocean)
print('{0} of {1} points in ocean'.format(ocean_samples, samples))

# Plot the detected sea points
map.drawcoastlines()
map.scatter(x.reshape(samples)[ocean.mask],
            y.reshape(samples)[ocean.mask],
            s=0.1, color='c', edgecolors='none')
plt.show()

The result I get:

80928 of 239888 points in ocean

Unfortunately, as you can see on the plot below, the resolution offered by maskoceans() is a bit coarse, even at the highest resolution and smallest grid size available.

Strictly speaking, if you are mainly interested in measuring area, it would be better to choose an equal-area projection, though it makes little difference at this zoom level.