Mostly so I myself can remember how to do it, here is how to create a random geotiff with GDAL in Python
Note: the width and height are given in opposite order in the GDAL raster and numpy arrays!
width = 4000
height = 3000
format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.Create( "test.tiff", width, height, 1, gdal.GDT_Byte )
dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
srs = osr.SpatialReference()
dst_ds.SetProjection( srs.ExportToWkt() )
raster = numpy.zeros( (height, width), dtype=numpy.uint32 )
color_range = 2**8
seed = math.pi**10
for i in range(height):
for j in range(width):
color = (seed*i*j) % color_range
raster[i][j] = color
dst_ds.GetRasterBand(1).WriteArray( raster )
It’s kind of slow, so perhaps the operation can be speeded up somehow? The result looks kind of nice though (image created with width and height both 4000):
So, not completely random.
pyproj is a Python interface to the PROJ.4 (http://trac.osgeo.org/proj/) functions. It allows you to transform coordinates between coordinate systems.
Install pyproj using pip
Assuming you’re using pip to install Python packages:
Here is a map application that illustrates very well the different tile schemes used (TMS, QuadTree, Google Maps):
Online documentation for GDAL/OGR Python is sparse. Here I show some recommended ways of learning more about GDAL/OGR in Python.
Using Python interpreter
You can learn about GDAL and OGR from inside the Python interpreter.
Start python interpreter:
Import the modules:
from osgeo import gdal,ogr,osr
Learn about the modules using ‘help’ and ‘dir’ in Python. These built-in functions work for any type of object (module, class, functions etc):
help(ogr) # display help for the ogr module
dir(ogr) # what's contained in ogr module?
help(ogr.Geometry) # display help for ogr.Geometry class
dir(ogr.Geometry) # show contents of an object, like functions on a class
gis.stackexchange.com is a Q&A site about GIS. Many people here know about GDAL and ORG, also about the Python bindings.
Some good questions I’ve found:
Transformation using OGR in Python
Will add more as I find them…
GDAL is not very well documented, and a complete Python API documentation site is not something I’ve found. So here are the best tutorials I’ve found so far for GDAL + Python.
GDAL Python samples on osgeo.org
Simple CSV file import
You have a CSV file called “data.csv”. It has a header line, and is delimited using “;”. You want to import it into Postgres and a table called “your_table”:
Create the database table. Set column-types so the string fields in the CSV file, can be cast to values in columns.
CREATE TABLE your_table
-- Your columns
Assuming a shapefile called myshapefile.shp, a table mytable in schema xyz, in a PostGIS enabled database called mydb on localhost. The table is owned by user dbuser who has password “secret”.
shp2pgsql myshapefile -I xyz.mytable > statements.sql
psql -d mydb -h localhost -U dbuser -f statements.sql
This tip and many more can be read in Making Maps Fast.
This is even easier with ogr2ogr:
ogr2ogr -f "PostgreSQL" PG:"host=localhost user=dbuser dbname=mydb password=secret" -lco SCHEMA=xyz myshapefile.shp
Google Fusion Tables is a service for storing tabular data online. It has an API that allows you to make an SQL-like request for data from an online table, like the following: sql=SELECT * FROM 1906375 (clicking the link will download a small CSV file)
In a previous post I tried generating OpenStreetMap tiles using GeoServer. It ended when I couldn’t find a good style (SLD) to apply the OSM layer.
The (new) idea, that also failed miserably
In this post I tried to use Mapnik to generating OSM map tiles in EPSG:25832. It failed mainly because the Python scripts published by OSM for generating tiles don’t support epsg:25832 out of the box. Mapnik is however an obvious choice for OSM because:
- Mapnik can read .osm files
- There is a comprehensive style for Mapnik, that is being maintained by OSM
Warning: This is a description of how to create a OpenStreetMap WMS with GeoServer. It works fine up to the point where the layers published as an unstyled WMS. This is where I’ve not been able to produce a good result, because of lack of a good Styled Layer Descriptor (SLD). If you have hints about a good SLD, feel welcome to submit a comment!
The idea is/was to first create a good general purpose OpenStreetMap WMS, and then use GeoWebCache to generate tiles from this WMS source in a custom projection, epsg:25832 in our case.