Tag Info

Hot answers tagged

17

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 ...


10

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 = ...


9

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" ...


8

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 <- ...


7

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 ...


5

For a complete history of the evolution from EPSG:900913 to EPSG 3785 and finally EPSG:3857, look at AnĂ¡lisis de Google Maps (in Spanish) and the conclusions are: The projection parameters defined for the EPSG: 900913 or 3785 describe a Mercator projection with geographic coordinates defined on a spherical model of the Earth, 6378.137m ratio, as specified ...


5

As the geometry is in the WTK-format and is of the type polygon, you will specify this in your vrt-file. So your vrt file should look someting like this: <OGRVRTDataSource> <OGRVRTLayer name="parcel"> <SrcDataSource>parcel.csv</SrcDataSource> <GeometryType>wkbPolygon</GeometryType> ...


5

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 ...


5

Avoiding the spatial query Since you noted buffering is computationally expensive and may be holding you back, consider this approach: Start looping through each point and round off your lat long point to a decimal place within your buffer (i.e. if your lat/long is 12.3456789/12.3456789 then get all points that begin with a lat/long of 12.34567/12.34567 or ...


5

You can solve with two chained VRT files and a bit of OGR SQL. The first VRT (e.g. remapped_csv.vrt) is: <OGRVRTDataSource> <OGRVRTLayer name="remapped_csv"> <SrcDataSource>test.csv</SrcDataSource> <SrcSQL>SELECT *, SUBSTR(latlon,2,5) AS lat, SUBSTR(latlon,9,12) AS lon FROM test</SrcSQL> ...


5

Within the database, geometries are stored on disk in a format only used by the PostGIS program. In order for external programs to insert and retrieve useful geometries, they need to be converted into a format that other applications can understand. Fortunately, PostGIS supports emitting and consuming geometries in a large number of formats: from ...


4

ESRI Shapefiles are more widely used, but have some limitations, including that the attribute names are limited in length. Shapefiles can be edited by the wonderful QGIS. MapInfo files are less widely used, and have a limitation with respect to something called bounds, where the X, Y coordinates are stored as an integer between two doubles. This limits the ...


4

Pyproj is a Python C extension which uses the PROJ4 library and osgeo.ogr is a Python C extension which uses the PROJ4 library. If you're only considering coordinate projection, in the fairest test they'd be almost equal. Pyproj's transform can operate on arrays of coordinate values, so you only need to call it once per line or ring instead of for every ...


4

It is a pure Python problem and not an GDAL/OGR problem. But if you want to program in Python,you need to understand the difference between an iterator and an iterable (many, many references on the Web, but one of the best is Loop Like A Native by Ned Batchelder). These are fundamental concepts in Python. "The important operation on an iterable is ...


4

The centroid function returns the centroid of the vertices within a geometry (see the documentation). It was originally only intended for use with polygons but can now be used with other geometries. So, the start-point of a line gives no indication of where the centroid of its vertices should be and is therefore irrelevant in your example above (I'm not ...


4

I guess you might want to append (>>) the projection to list.txt and not assign (>), otherwise list.txt will be overwritten at each iteration of your for loop. For /R %f in (*.tab) do ogrinfo -al -so %~dpnxf >> list.txt


4

I've gotten this working on my 10.1 system, at least somewhat (I can run import ogr from the ArcGIS python window and it appears to be working). First I attempted to get my system as clean as possible. I uninstalled previous versions of GDAL (from OSGeo4W). I do still have PostgreSQL+PostGIS which has it's own version of GDAL, but I can't uninstall that (I ...


4

As it already has been mentioned in comments to your questions, MapInfo Professional does require the table to be added to the MapInfo.MapInfo_MapCatalog. This table tells MapInfo Professional how to handle and display the spatial part of the table. You do however not need to use EasyLoader to load your table - you can do it with other tools, too. But you ...


4

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 ...


4

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 ...


3

It appears your geometry is unknown to ogr. It doesn't see your geometries as one type and reports: kyle@kyle-workstation:Downloads$ ogrinfo -so -al Example.dxf INFO: Open of `Example.dxf' using driver `DXF' successful. Layer name: entities Geometry: Unknown (any) Feature Count: 7886 Extent: (2481827.566011, 337655.197584) - (3104773.310852, ...


3

In general GDAL/OGR don't throw errors in Python when they can't retrieve a feature object. Instead they return None objects. This means that you have to manually test for a feature existing when you retrieve it. Secondly, if you're using GetNextFeature, OGR reads only in one direction, like if you were reading a text file, so to go back to the beginning of ...


3

If using GDAL 1.10+ the python bindings allow you to specify a python callable as an error handler. However, these error handlers appear to be called in a separate thread and any exceptions raised do not propagate back to the main thread. So something like this won't work: from osgeo import gdal def error_handler(err_level, err_no, err_msg): if ...


3

Taking the information from the above, here's how you would actually get the ogr2ogr command working like it should add this to your ~/.bash_profile script: export PATH=/Library/Frameworks/GDAL.framework/Versions/1.[YOURVERSIONHERE]/Programs/:$PATH the above code also adds a bunch of other gdal libraries to your path as well


3

Early in your program, you've done this: for feature1 in feat1Layer: geomfeat1 = feature1.GetGeometryRef() After which the feat1Layer "cursor" is at the end of the file. Later, when you call inFeature = feat1Layer.GetNextFeature() GetNextFeature() will return None. And if it wasn't None, the program would fail at the end of the loop with a ...


3

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 ...


3

I tried this example and I got different results. For a shapefile containing a layer with 246 features, the running time was 304usec per loop on GetFeatureCount(), while it was 315usec per loop on len(layer). Timeit used 1000 loops to evaluate this. This was done using: testFeatureCount.py import osgeo.ogr def testFeature(): shapefile = ...


3

len() is a pure Python function that return the length (the number of items) of an object (whatever the object) -> in this case the number of features. .GetFeatureCount() did the same thing (number of features in the layer) but has some additional functions compared to the simple len(). It is based on the CPP OGRLayer::GetFeatureCount() function (in the ...


3

ogrinfo as nothing to do with Visual Studio. It is a command line tool (Windows, Linux, Mac OS X) So you use it in the Windows console or in the Linux and Mac OS X shell in the same way. ogrinfo yourshape.shp with a lot of options: (ogrinfo) And, there are many, many tutorials in Gis.StackExchange : search ogrinfo and on the Web....


3

Quoting the shapefile specification (page 4): All the non-Null shapes in a shapefile are required to be of the same shape type. You may choose point OR multipoint OR polygon, but you cannot use multiple geometry types in a single shapefile. Multi-part geometries are permitted in all flavors of multipoint, polyline, and polygon shapefiles.



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