Skip to content

OGR to reproject, modify Shapefiles

nvkelso edited this page Apr 1, 2013 · 31 revisions

Install

GDAL includes OGR.

Ubuntu 11.10 - Source

Check to see if you have binutils: dpkg -s binutils

If you don't have binutils, sudo apt-get install binutils

sudo apt-get install gdal-bin

sudo apt-get install libproj-dev

Upgrading to new versions »

Mac OS X

Install the KyngChaos GDAL frameworks for your version of Mac OS X.

OGR

Usage

VPRO reprojection:

to geographic

ogr2ogr -t_srs EPSG:4326 outputwith4236.shp input.shp

ogr2ogr -s_srs "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs" -t_srs EPSG:4326 output_shapefile.shp input_shapefile.shp

to mercator:

OGR2OGR -f "ESRI Shapefile" -s_srs "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs" -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over" /Users/nvkelso/github/golden-ratio/scripts/900913/flickr_bbox_merc.shp /Users/nvkelso/github/golden-ratio/scripts/flickr_bbox_join8_redux_poly_area1.shp


ogr2ogr -s_srs "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs" -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"  output_shapefile.shp input_shapefile.shp

Rename SHP columns:

ogrinfo -so z7-labels.shp z7-labels

ogr2ogr -select 'name, admin1 cod as admin1_cod, point size as point_size , type_rank, name, geonameid, font file as font_file, zoom as zoom_dymo, country co as country_co, font size as font_size, capital, asciiname, type, population’

Sort the OGR columns:

Important for drawing features in the right order (especially town labels, also important for how roads stack).

ogr2ogr -sql "SELECT FEATURECLA, NAME, AGGLOMERA2, BILINGUAL, POP_MAX, POPRANK, LABELRANK, SCALERANK, ADM0_A3, GEONAMEID FROM pop12_ru_merc3 ORDER BY POP_MAX DESC" shp/pop12_ru_merc3{-sorted,}.shp


ogr2ogr -sql "SELECT * FROM dma_cities_a_b_c_d_join ORDER BY POP_MAX DESC" dma_cities_a_b_c_d_join{_sorted,}.shp

ogr2ogr -sql "SELECT * FROM airports_merc ORDER BY natlScale DESC" airports_merc{_sorted,}.shp

Stats on your values

Use -sql flag for extra power. Docs

Shows me I have a couple hundred country codes and what they are:

ogrinfo -sql "select COUNT(DISTINCT ADM0_A3) from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple

ogrinfo -sql "select DISTINCT ADM0_A3 from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple

Shows me there are just a handful of scale ranks:

ogrinfo -sql "select COUNT(DISTINCT SCALERANK) from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple

Describes the actual values in that column:

ogrinfo -sql "select DISTINCT SCALERANK from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple

Find out the min, max, avg values in a number column:

ogrinfo -sql "select MIN(SCALERANK) from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple
ogrinfo -sql "select MAX(SCALERANK) from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple
ogrinfo -sql "select AVG(SCALERANK) from ne_10m_populated_places_simple"  ne_10m_populated_places_simple.shp ne_10m_populated_places_simple

Exclude high or low extreme values:

SELECT * from polylayer WHERE prop_value > 220000.0 ORDER BY prop_value DESC

Densify your linework:

Note: Uses the units of your data store, here 0.5 degrees assuming Geographic.

ogr2ogr -segmentize 0.5 <out_file> <in_file>

Clipping linework with a shape:

ogr2ogr -clipsrc clipping_polygon.shp output.shp input.shp

Clipping linework with a boundingbox:

tk tk tk

OGR in Python

Has got to be the worst documented API ever.

ChrisG at USU.edu has a good tutorial

Want the power of OGR in Python? Try this other tools: Fiona from @sgillies wrapping OGR and the awesome shapely.

Gotcha's

Ogr2Ogr turns KML into Polyline MZ files instead of simple Polylines. In ArcMap, you'll need to convert it down to do topology and other analysis operations: http://forums.esri.com/Thread.asp?c=93&f=993&t=228662

GDAL

http://trac.osgeo.org/gdal/wiki/FAQRaster

How to make GeoTIFF from non-geospatial raster?

One option is to use gdal_translate to assign the geo-referenced bounds and generate .tfw file:

gdal_translate -a_ullr ulx uly lrx lry src_dataset dst_dataset

or with SRS assignment:

gdal_translate -a_srs EPSG:4326 -a_ullr ulx uly lrx lry src_dataset dst_dataset

How do I use gdal_translate to extract or clip a sub-section of a raster?

Gdal_translate is designed to convert to and from a variety of raster formats, but it can also perform helpful geoprocessing operations during conversion.

If you would like to extract a sub-section of a raster you can use the -srcwin or -projwin options. In gdal terminology these are "subsetting" operations that allow you to select a "subwindow" to copy from the source dataset into the destination dataset.

Here is an example of using gdal_translate on NAIP orthophotography in sid format to select a small subwindow that shows Blakely Island, WA:

gdal_translate -projwin 510286 5385025 518708 5373405 ortho_1-1_1n_s_wa055_2006_1.sid naip_ortho_blakely_island.tif

This example uses the -projwin option which accepts the bounding coordinates in the projected coordinates rather than in pixels (-srcwin). Gdal_translate -projwin needs the upper left x coordinate, the upper left y coordinate, the lower right x coordinate, and the lower right y coordinate. The naip imagery in this example is in NAD 83 Utm 10, so to get these bounding coordinates I simply loaded up the index shapefile that comes packaged with naip imagery in Quantum GIS and read the screen coordinates to form my extent.

Another way

http://www.gdal.org/gdalwarp.html

gdalwarp -s_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -t_srs EPSG:4326 -te -123.04962 37.54675 -121.81091 38.23387 -r near raw_data/usgs_nlcd/nlcd2006_landcover_4-20-11_se5.img nlcd2006_mill_valley.tif

Note: Currently, clipping raster images using vector extent polygons is not supported but is under discussion (see http://trac.osgeo.org/gdal/ticket/1599). However it is fairly easy to get the extents of a given shapefile and convert those coordinate pairs into the format needed by gdal_translate without manually reading the extents from another application like qgis. Say you have a shapefile named clipping_mask.shp use ogrinfo to get the extents:

Note the use of the pipe and grep command is optional(| grep Extent), but is a slick way to limit the info reported by ogrinfo to just what you need in this case

ogrinfo clipping_mask.shp -so -al | grep Extent

# which gives the extent as xMin,yMin, xMax, yMax:
Extent: (268596, 5362330) - (278396, 5376592)
# which is (xMin,yMin) - (xMax,yMax)

Then copy and paste that text to create your gdal_translate clipping command:

# -projwin's ulx uly lrx lry is equivalent to xMin, yMax, xMax, yMin so just switch the Y coordinates
# For the above Extent that would turn into:

gdal_translate -projwin 268596 5376592 278396 5362330 src_dataset dst_dataset

How to create or modify an image color table ?

The easiest non-scripting way to create or modify a color palette is to translate your file to VRT (XML) format, and edit the color table in the XML file.

# create the vrt
gdal_translate -of VRT your.tif your.vrt

## --- edit colour table here ---

# create final fixed image
gdal_translate your.vrt your_fixedup.tif

The .vrt file doesn't actually include the image data, it just references back to your.tif (or other supported format). This general approach can be used to modify metadata, color tables, gcps and other things for which there is no easy to do set them with gdal_translate. There are other mechanisms to do this via scripting or c/c++ programming.

Here is an example for USGS DRGs, see Virtual Format Tutorial for detailed syntax options. The entries are in palette order and correspond to red, green, blue & alpha channels.

<ColorInterp>Palette</ColorInterp>
<ColorTable>
    <Entry c1="0" c2="0" c3="0" c4="255"/>
    <Entry c1="255" c2="255" c3="255" c4="255"/>
    <Entry c1="0" c2="151" c3="164" c4="255"/>
    <Entry c1="203" c2="0" c3="23" c4="255"/>
    <Entry c1="131" c2="66" c3="37" c4="255"/>
    <Entry c1="201" c2="234" c3="157" c4="255"/>
    <Entry c1="137" c2="51" c3="128" c4="255"/>
    <Entry c1="255" c2="234" c3="0" c4="255"/>
    <Entry c1="167" c2="226" c3="226" c4="255"/>
    <Entry c1="255" c2="184" c3="184" c4="255"/>
    <Entry c1="218" c2="179" c3="214" c4="255"/>
    <Entry c1="209" c2="209" c3="209" c4="255"/>
    <Entry c1="207" c2="164" c3="142" c4="255"/>
</ColorTable>

Import into PostGIS directly, using Python bindings

See this Gist from @migurski: https://gist.github.com/3778300

Use Python ORG bindings for direct import into PostGIS

http://www.gdal.org/ogr/drv_pg.html

Data joins

More details »

OGR SQL supports a limited form of one to one JOIN. This allows records from a secondary table to be looked up based on a shared key between it and the primary table being queried. For instance, a table of city locations might include a nation_id column that can be used as a reference into a secondary nation table to fetch a nation name. A joined query might look like:

SELECT city.*, nation.name FROM city 
     LEFT JOIN nation ON city.nation_id = nation.id
Something went wrong with that request. Please try again.