-
Notifications
You must be signed in to change notification settings - Fork 50
OGR to reproject, modify Shapefiles
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
Mac OS X
Install the KyngChaos GDAL frameworks for your version of Mac OS X.
#OGR
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
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
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.
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
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