attribute: Phillie Casablanca

Category: GIS (1)

Understanding raster, basic GIS concepts and the python gdal library



I've been learning about how to handle *raster* files and getting comfortable with the python libraries that are available to manipulate them.
(I'm still learning, but hopefully this will help someone like me who is just getting started)

First thing you need to know about common *decimal degree* latitude/longitude notation is that it is **NOT** a linear representation, but *angular* (defining a specific point on a sphere-like shape). So performing `cartisian coordinate <http://en.wikipedia.org/wiki/Coordinate_ssytem>`_ (where the distance from (1,1) to (1,2) and the distance from (1,1) to (1,0) is an equal distance) operations is going to give you very accurate results.

Next, the common *decimal degree* latitude/longitude values you see use a specific **geographic coordinate system** (Commonly *WGS84*). In addition to **geographic**, **projected coordinate systems** are also used.

This quote from the `OGR Projections Tutorial <http://www.gdal.org/ogr/osr_tutorial.html>`_ helped me to understand what's going on here:

>> There are two primary kinds of coordinate systems. The first is **geographic** (positions are measured in *long/lat*) and the second is **projected** (such as UTM - positions are measured in *meters* or *feet*).

>> A **projected** coordinate system (such as UTM, Lambert Conformal Conic, etc) requires an underlying **geographic** coordinate system as well as a definition for the *projection transform* used to translate between linear positions (in meters or feet) and angular long/lat positions.

In my examples here I'll be using **UTM** for the *projected* coordinate system and **WGS84** for the *geographic*. I'll be calling the "projection transform", *geo_transform* in my code below.

.. note::

UTM is split into multiple *zones* and may be difficult to work with in places where you will be spanning zones. The "Google" projection, *900913*, may be better for simple manipulations.

First some simple coordinate conversion. When doing manipulation of data you're going to need to do some conversions, here's how to do it with the python `GDAL <http://pypi.python.org/pypi/GDAL/>`_ package:

import osr

def transform_utm_to_wgs84(easting, northing, zone):
utm_coordinate_system = osr.SpatialReference()
utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
is_northern = northing > 0
utm_coordinate_system.SetUTM(zone, is_northern)

wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system

# create transform component
utm_to_wgs84_geo_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (, )
return utm_to_wgs84_geo_transform.TransformPoint(easting, northing, 0) # returns lon, lat, altitude

def transform_wgs84_to_utm(lon, lat):
def get_utm_zone(longitude):
return (int(1+(longitude+180.0)/6.0))


def is_northern(latitude):
"""
Determines if given latitude is a northern for UTM
"""
if (latitude < 0.0):
return 0
else:
return 1

utm_coordinate_system = osr.SpatialReference()
utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
utm_coordinate_system.SetUTM(get_utm_zone(lon), is_northern(lat))

wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system

# create transform component
wgs84_to_utm_geo_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (, )
return wgs84_to_utm_geo_transform.TransformPoint(lon, lat, 0) # returns easting, northing, altitude

As you may have noticed, *easting* and *northing* are the *UTM* terms for *x*, and *y*.

Now, why would you want to use *rasters*? Well, rasters are used as a way to aggregate data in a given area. I've commonly called this "*binning*", but in raster terms the "*bin*" elements of a raster are refered to as a "**cell**" or "**pixel**". The *raster* structure/object model in gdal is like this:

- A *raster* consists of a *dataset*.

- A *dataset* consists of one or more *bands*.

- A *band* contains multiple *cells*/*pixels*/*bins*.

- A *cell*/*pixel*/*bin* is a location with a *value*.

Now for the code. I think there may be a good way to create a *RasterBuilder* class to handle the raster functions, but I'll leave that to you.

/

import numpy
from osgeo import gdal
from osgeo import osr
gdal.AllRegister() # register all drivers
gdal.UseExceptions()
from shapely.geometry import MultiPoint # used to get extend/bounds

def get_point_bounds(points):
mp_points = MultiPoint(points) #from shapely.geometry import MultiPoint
(min_lon, min_lat, max_lon, max_lat) = mp_points.bounds
return min_lon, min_lat, max_lon, max_lat

def get_raster_size(minx, miny, maxx, maxy, cell_width, cell_height):
"""
Determine the number of rows/columns given the bounds of the point data and the desired cell size
"""
cols = int((maxx - minx) / cell_width)
rows = int((maxy - miny) / abs(cell_height))
return cols, rows

def lonlat_to_pixel(lon, lat, inverse_geo_transform):
"""
Translates the given lon, lat to the grid pixel coordinates in data array (zero start)
"""
# transform to utm
utm_x, utm_y, utm_z = transform_wgs84_to_utm(lon, lat)

# apply inverse geo tranform
pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y)
pixel_x = int(pixel_x) - 1 # adjust to 0 start for array
pixel_y = int(pixel_y) - 1 # adjust to 0 start for array

return pixel_x, abs(pixel_y) # y pixel is likly a negative value given geo_transform

def create_raster(points, filename="test.tiff"):
"""
points == { (lon, lat): value, ...}
lon/lat values in WGS84
"""
# create empty raster
driver = gdal.GetDriverByName(OUTPUT_FORMAT)
number_of_bands = 1
band_type = gdal.GDT_Float32

x_rotation = 0
y_rotation = 0
cell_width_meters = 50.0
cell_height_meters = 50.0

(min_lon, min_lat, max_lon, max_lat) = get_point_bounds(points) # retrieve bounds for point data

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
srs.SetUTM( get_utm_zone(min_lon), is_northern(max_lat)) # Set projected coordinate system to handle meters

# create transforms for point conversion
wgs84_coordinate_system = srs.CloneGeogCS() # clone only the geographic coordinate system


wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, srs)
utm_to_wgs84_transform = osr.CoordinateTransformation(srs, wgs84_coordinate_system)

# convert to UTM
top_left_x, top_left_y, z = transform_wgs84_to_utm(min_lon, max_lat)
lower_right_x, lower_right_y, z = transform_wgs84_to_utm(max_lon, min_lat)

cols, rows = get_raster_size(top_left_x,
lower_right_y,
lower_right_x,
top_left_y,
cell_width_meters,
cell_height_meters)
dataset = driver.Create(filename, cols, rows, number_of_bands, band_type)

# GeoTransform parameters
# --> need to know the area that will be covered to define the geo tranform
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, -cell_height_meters ] # cell height must be negative (-) to apply image space to map
dataset.SetGeoTransform(geo_transform)
dataset.SetProjection(srs.ExportToWkt())
inverse_geo_transform = gdal.InvGeoTransform(self.geo_transform)[1] # for mapping lat/lon to pixel

# get the empty raster data array
band = dataset.GetRasterBand(1) # 1 == band index value
data = band.ReadAsArray(0, 0, cols, rows).astype(numpy.cfloat)

# populate array values for output
for lon, lat in points.keys():
value = points[(lon, lat)]
# apply value to array
pixel_x, pixel_y = lonlat_to_pixel(lon, lat, inverse_geo_transform)
data[pixel_x][pixel_y] = value


# write the updated data array to file
band.WriteArray(data, 0, 0)
band.SetNoDataValue(NULL_VALUE)
band.FlushCache()

# set dataset to None to "close" file
dataset = None


monkut // May 2, 2012 // 6:15 p.m.