# How to map all 0.001° graticules that contain land?

I'm trying to map all 0.001° graticules (across the planet) that contain land.

My first attempt naively tried to brute force this by literally looping over the entire planet in 0.001° increments and checking if each point contained land or not. Should've done the math first -- not optimal :(.

My second attempt tried narrow down the scope by first finding the bounding boxes for known land areas and then doing the same "loop and check" technique over that. I used the "land polygons" from NaturalEarth, however, the bounding boxes I was able to generate from those polygons weren't granular enough.

I'm thinking next I'll try to find a larger set (of smaller polygons) so the bounding boxes are more "accurate"... maybe countries instead of land masses?

Either case, I'm sure this approach will eventually work, but, I'm thinking there's gotta be something better/more optimal. (Also, it should probably be noted that I'm incredibly new to GIS related things :).

Thoughts? Is there a way to do this via a query inside of PostGIS (which is what I'm currently using to store the data)? Or is there a better alternative?

Thanks, and let me know if you need any extra clarification.

• By "graticule" would you perhaps mean a point whose (lat, lon) coordinates are whole multiples of 0.001 degree? By "point contained land" do you mean that such a point is within a land area? Are you aware there are approximately 20 billion (2E+10) such points? And that they are not equally spaced? These questions raise so many issues of practicability that perhaps it would be more constructive to back up a little and explain why you are trying to do this: perhaps there is a much more efficient way to achieve your objectives. – whuber Aug 5 '14 at 19:27
• Another thing to consider here is that whatever the source of the world map that you're basing the analysis on is, it's likely at a much coarser resolution (scale) than the scale that you're trying to extract information from it. Doing so is usually a bad idea unless it is well justified. – WhiteboxDev Aug 5 '14 at 20:53

Burning polygons into a raster is by far the most efficient way to classify a graticule as containing land (or not). The `ALL_TOUCHED=TRUE` rasterizing option ensures that a graticule location is "contained" if even a small part of a polygon touches it. The processing should take a few minutes to do.

``````from osgeo import ogr, gdal

# Shapefile, PG: connection, or whatever, with polygons of land/not-land
ogr_ds = ogr.Open(datasource_to_vector_data)
ly = ogr_ds.GetLayer()

# Create a 1-bit GeoTIFF raster, which is a grid of land/not-land
drv = gdal.GetDriverByName('GTiff')
ds = drv.Create('earth.tif', 360000, 180000, 1, gdal.GDT_Byte, ['NBITS=1'])
ds.SetGeoTransform([-180, 0.001, 0.0, 90.0, 0.0, -0.001])
band = ds.GetRasterBand(1)
band.Fill(0)  # 0 = not-land; burn 1 = land
gdal.RasterizeLayer(ds, [1], ly, None, None, [1], ['ALL_TOUCHED=TRUE'])
ly = ds = band = ogr_ds = None  # save, close
``````

The resulting file is around 8 GB, and can be used to determine if each location is land (1) or not (0).

The raster result can be loaded into PostGIS using the `1BB` 1-bit boolean pixel type. Further note that all of these steps can be done within PostGIS using `ST_AsRaster`, which uses the same GDAL engine for rasterizing and burning.

• I was just about to add a 360000x180000 1-bit raster (8Gb, uncompressed) to my answer. It's still a fearsome pixel count (and loses precision within individual tiles), but you can use any one of a number of high-resolution vector sources for population. – Vince Aug 6 '14 at 1:11
• What's the count of '1' pixels? – Vince Aug 6 '14 at 13:37
• @Vince That count may be practically irrelevant. Internally the raster is likely stored as rows which can be extremely efficiently encoded with a lossless run-length encoding. Assuming an average of, say, 50 changes of pixel across each row, the total file might need only about (180/0.001)*50 = 9 million (long) integers to represent all the information with no loss. Furthermore, a "smart" GIS could perform many raster calculations without having to expand this compressed representation into an uncompressed form. – whuber Aug 6 '14 at 14:29
• I was just curious to see if we ball-parked the count correctly. – Vince Aug 6 '14 at 15:23
• @Vince reading the raster band to a Numpy array takes just over 60 GB of RAM and slightly more time to process. I counted 21461512201 or 33.1% of total as land. Mind you, my vector datasource has a large-scale resolution, and isn't fantastic for this resolution. – Mike T Aug 6 '14 at 22:54

So you realize you're trying to identify roughly 23 of 64 billion tiles? Assuming that each query takes 10 milliseconds (which is still too low), you're looking at 20+ CPU years. Since large swaths of the Earth are all land or all water, you can cut the corner on these, but assuming this reduces the problem to only 1/10th of 64,800 one-degree grids, you've still likely got 2 CPU-years of effort just to identify which 23 billion tiles are to be mapped. Constructing 23B maps for each tile with land is probably 1-3 seconds of effort to generate a 4k image, and you're back into seven hundred to two thousand CPU-years to generate 86 petabytes of images.

If you use memory-based techniques, with a mean query time of 50 microseconds, you could generate the tile list in only 38 CPU-days, and if all you're mapping is land and water, it might only take 100 milliseconds each to generate the 6.5 billion 4k images which have both land and water, which would still bring you to 70 CPU-years and 24 petabytes.

My thought is that this is not a reasonable project with current technology.