Hot answers tagged

35

You can use the ogr2ogr utility which is packaged with the gdal command line tools. Use the -sql option as follows: ogr2ogr outputfile.shp inputfile.shp -sql "SELECT oldfield1 AS newfield1, oldfield2 AS newfield2 from inputfile" As an added bonus, you can convert the data into a different format at the same time, or filter your data by specifying a where ...


32

You miss a minus sign before where and the select is not necessary, so it should be: ogr2ogr -where ID="1" outfile.shp infile.shp or if you have to do more complex query on your input data: ogr2ogr -sql "SELECT * FROM infile WHERE ID='1'" outfile.shp infile.shp If ID is a field of Integer type, substitute ID='1' with ID=1. Notes: -f "ESRI Shapefile" is ...


31

OpenGIS Simple Features Reference Implementation. You can read a short statement about the history of the abbreviation among some other GDAL related facts in the GDAL FAQ. The relevant paragraph states: OGR used to stand for OpenGIS Simple Features Reference Implementation. However, since OGR is not fully compliant with the OpenGIS Simple Feature ...


30

In pure Python, without using the subprocess module (os.system is deprecated) to call ogr2ogr or shp2pgsql, for example): 1) with ogr Append a shapefile to a postgis table using the GDAL/OGR Python interface 2) with ogr and psycopg2 from the book Python Geospatial Development (Eric Westra), Chapter 7, p.219 import os.path import psycopg2 import ...


29

ogr2ogr is part of the Geospatial Data Abstraction Library (GDAL). Get homebrew from http://brew.sh brew install gdal


29

Shapefiles have no type MultiPolygon (type = Polygon), but they support them anyway (all rings are stored in one feature = list of polygons, look at Converting huge multipolygon to polygons) The problem If I open a MultiPolygon shapefile, the geometry is 'Polygon' multipolys = fiona.open("multipol.shp") multipolys.schema {'geometry': 'Polygon', '...


27

Only adding this because I tried using the kyng chaos tools, but on my Mac OS X machine I was able to very, very easily install this with Anaconda conda install gdal Posting in case anyone finds this again - I realize the original post is 3 years old.


26

The second argument in Open specifies if the data can be updated (written to). Try: dataSource = driver.Open(fn,1)


24

The ogr2ogr utility supports a limited sql syntax. You can join your CSV to the shapefile using something like the following: ogr2ogr -sql "select inshape.*, joincsv.* from inshape left join 'joincsv.csv'.joincsv on inshape.GISJOIN = joincsv.GISJOIN" shape_join.shp inshape.shp


20

As the name says, a SpatialPolygonsDataFrame is basically just a SpatialPolygons object with data attached (the attribute table). The data must have at least as many rows as there are features library(rgdal) ob <- SpatialPolygons(..)# Your SpatialPolygons Object spp <- SpatialPolygonsDataFrame(ob,data=as.data.frame("yourData"),proj4string=CRS("+...


19

ogrinfo can shorten the output considerably using the -so flag. -so: Summary Only: supress listing of features, show only the summary information like projection, schema, feature count and extents. So ogrinfo -ro -so file.shp should give a summary of the metadata. And -al: List all features of all layers (used instead of having to give layer names ...


19

This tutorial shows how to import geopackges into postgres with ogr2ogr like: ogr2ogr -f PostgreSQL "PG:dbname=gadm" DNK_adm.gpkg Note importing into PostgreSQL like this will convert the case of the field names in your GeoPackage data to lowercase in the PostgreSQL table. To maintain case you would need to use the lco option like: ogr2ogr -f PostgreSQL "...


18

The following script determines the bounding box of a raster and creates based on the bounding box a geometry. import ogr, gdal raster = gdal.Open('sample.tif') vector = ogr.Open('sample.shp') # Get raster geometry transform = raster.GetGeoTransform() pixelWidth = transform[1] pixelHeight = transform[5] cols = raster.RasterXSize rows = raster.RasterYSize ...


18

There's a special field in OGR SQL called OGR_GEOM_AREA which returns the area of the feature's geometry: ogrinfo -sql "SELECT SUM(OGR_GEOM_AREA) AS TOTAL_AREA FROM myshapefile" myshapefile.shp where TOTAL_AREA unit of measure depends by the layer SRS (read the comments below).


17

I had to solve the same problem today, so here is my answer, which gives a complete solution. I have a lineWKT.csv file stored in F:\Data\ folder, with the data like this: id,gm 0,"LINESTRING (30 10 0, 10 30 0, 40 40 5)" I have a test.vrt file like this: <OGRVRTDataSource> <OGRVRTLayer name="lineWKT"> <SrcDataSource>F:\Data\...


17

So far brew install gdal --HEAD resulted in $ ogr2ogr --version GDAL 2.1.0dev, released 2015/99/99


14

Using GDAL >= 1.10.0 compiled with SQLite and SpatiaLite: ogr2ogr data_shifted.shp data.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,1,10) FROM data" where shiftX = 1 and shiftY = 10.


13

There are some good libraries out there that do most of the heavy lifting for you. Example using [shapely][1] in python. import random from shapely.geometry import Polygon, Point def get_random_point_in_polygon(poly): minx, miny, maxx, maxy = poly.bounds while True: p = Point(random.uniform(minx, maxx), random.uniform(miny, maxy)) ...


13

Because UPDATE is not supported in OGR SQL, as you stated in a comment, you should update the table using the SQLite SQL dialect available in GDAL >= 1.10 with SQLite and SpatiaLite support: ogrinfo $myfile -sql "ALTER TABLE $name ADD COLUMN code_num integer(3)" ogrinfo $myfile -dialect SQLite -sql "UPDATE $name SET code_num = CAST(code_06 AS integer(3))"


13

The best all round tool here is a raster calculator. gdal_calc is a GDAL raster calculator implemented in Python here, with some examples here. If you e.g. wants to keep values above +50: gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>50)" --NoDataValue=0 You can specify several files -A to -Z, where each of them get a corresponding ...


13

Looks like ogr2ogr has support for exactly this problem. The below is copied directly from that page: How do I flip coordinates when they are not in the expected order The EPSG has a recommanded order for geographic SRS where the coordinates tuples of a geometry must appear in the (latitude, longitude) order, whereas most GIS will properly display such ...


12

In order to get the coordinates in decimal degrees, the data needs to be reprojected to WGS84. import ogr, osr driver = ogr.GetDriverByName('ESRI Shapefile') shp = driver.Open('testpoint.shp', 0) lyr = shp.GetLayer() feat = lyr.GetNextFeature() geom = feat.GetGeometryRef() # Transform from Web Mercator to WGS84 sourceSR = lyr.GetSpatialRef() targetSR = ...


12

Use these commands: sudo apt-get install python-gdal sudo apt-get install gdal-bin


12

It is easy: from osgeo import ogr import os driver = ogr.GetDriverByName("ESRI Shapefile") if os.path.exists('your.shp'): driver.DeleteDataSource('your.shp')


12

Internally, PostGIS stores geometries in a binary specification, but it is queried and viewed outside as a hex-encoded string. There are two popular variations of well-known binary (WKB): EWKB (via ST_AsEWKB) - an extended WKB specification designed by PostGIS. OGC WKB (via ST_AsBinary) - specified by the OGC and ISO. For a while it was 2D-only, but later ...


12

The accepted answer is really useful, but I found that it was slow with a large-ish database. I believe it also limits your options when joining the data. My method now is to pull everything into SQLite (using a combination of csvkit and ogr2ogr): csvsql --db sqlite:///myjoindb.db --insert myjoincsv.csv ogr2ogr -append -f "SQLite" myjoindb.db myjoinshp.shp ...


12

At the following link all SQL functions that test spatial relationships are described: http://www.gaia-gis.it/spatialite-2.4.0/spatialite-sql-2.4.html#p12. The following spatial relationships can be used within a query ST_Equals - ST_Disjoint - ST_Touches - ST_Within - ST_Overlaps - ST_Crosses - ST_Intersects - ST_Contains - ST_Relate For this type of ...


10

Problems: 1: the outcome of UnionSpatialPolygons is a spatial polygon 2: converting the result back into a spatial polygon data frame is a real pain -a. you need a very exact data frame to attach to a spatial polygon -b. data you used for UnionSpatialPolygons has more rows than the output and is not formatted in the way that is needed. My (ugly) ...


10

Why not work globally ? calculate the distances between all points union the resulting lines pointx - pointy with a distance < 14m I will use Shapely, much easier for resolving these kinds of problems. You must iterate through all pairs of points to calculate the distance once (as distance point1-point2 = distance point2-point1). There are many ...


10

The Ogr function GetEnvelope() returns "a tuple (minX, maxX, minY, maxY)" (from here), but what you want (from what I can understand) is a Polygon describing the envelope/bbox? This is actually rather simple, as the tuple (minX, maxX, minY, maxY) is all you need to create a Polygon. Just create a Polygon based these, like so: from osgeo import ogr def ...


Only top voted, non community-wiki answers of a minimum length are eligible