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.