Tag Info

Hot answers tagged

5

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 vrt file like this: <OGRVRTDataSource> <OGRVRTLayer name="lineWKT"> ...


5

To process large files, you need to use a generator which only read/write one line at a time and Python has a command that does that: with. The with statement handles opening and closing the file, including if an exception is raised in the inner block. The for line in f treats the file object f as an iterable, which automatically uses buffered IO and ...


4

This is a GDAL/OGR Python gotcha. So when you use print inputlyr[0].GetGeometryRef(), the reference to the indexed feature has been cleaned up before geometry functions can be used. This is expected behaviour, and there is no roadmap to change it. The issue is that the GDAL/OGR Python bindings are over a decade old, and are tightly coupled with the ...


4

Some remarks because you mix many things ( if you want yo use your script outside the console of QGIS, install the Python module GDAL (osgeo) in your Python installation: it is installed in the Python version of QGIS if you are on Windows). 1) There are no Shapely geometries in your script, only PyQGIS geometries QgsPoint(x_min, y_min), ... ogr geometries ...


4

Python does not contain a separate library for GDAL but rather bindings to access the GDAL libraries. As GDAL is updated so are the affected bindings, ensuring that you'll have full access to GDALs functionality from Python. If you upgrade GDAL you will also upgrade the Python bindings (if the version you upgrade to supports Python bindings). If you build ...


3

Resolved! It's a module PATH Python console problem. In QGIS Brighton/apps/Python27/Lib/site-packages create osgeo.pth (Administrators) file and put # .pth file for the osgeo extensions osgeo start QGIS and install flowmapper ;) Clarifying, after create osgeo.pth file, open him and just add in first line the word osgeo, save and run QGIS! Done.


3

You're trying to use an OGR (vector) driver with GDAL (raster) tools. Here's a few lines of my working code that may help: char* BasePath = new char[FullPathMax]; // this does have a value before it's used OGRRegisterAll(); OGRDataSource *hDS = NULL; OGRSFDriver *Driver = NULL; hDS = OGRSFDriverRegistrar::Open(BasePath,FALSE,&Driver); As you can ...


3

Numbers in any programming language does not have "absolute" precision. As they must be represented as bits in the computer hardware its precision is limited. Anyway review this links, because maybe they can help you: http://comments.gmane.org/gmane.comp.gis.gdal.devel/19331 http://lists.osgeo.org/pipermail/gdal-dev/2011-August/029793.html


3

I ran into this exact problem because I hadn't extracted the zip file into the working directory. I was rushing to get the files converted and when I checked to make sure they were there, the "preview" that I got when I double-clicked on the zipped file tricked me into thinking they were available for conversion. Even if this wasn't the solution to your ...


3

ptrv/gpx2spatialite does this remarkably well, saving timestamps for all points and deriving speed and length data for tracks. It also won't import duplicate tracks, so you can feed it a huge pile for GPX files and it will munge them appropriately. Update: usage examples, as requested: Initialize new database: gpx2spatialite_create_db db.sqlite Add a ...


3

I used the ogr java bindings: public static void unitOfCRS(){ SpatialReference poSourceSRS = new SpatialReference(); // output: metre poSourceSRS.ImportFromEPSG(3068); System.out.println(poSourceSRS.GetAttrValue("UNIT")); // output: degree poSourceSRS.ImportFromEPSG(4326); ...


3

A python solution is fairly simple using ogr.GetDriverByName(in_format).Open(in_file). import sys from osgeo import ogr def main(in_file, in_format, out_file, out_format): in_ds = ogr.GetDriverByName(in_format).Open(in_file) out_ds = ogr.GetDriverByName(out_format).CopyDataSource(in_ds, out_file) if __name__ == '__main__': main(*sys.argv[1:]) ...


3

Assuming you have an older version of GDAL/OGR, you can use OGR SQL to cast the geom field to a geometry. For example: ogr2ogr -f "ESRI Shapefile" "sample.shp" "geo.csv" -sql "SELECT *, CAST(geom as geometry) FROM geo" This will create a new shapefile using the WKT data as the geometry. You can use the same query with ogrinfo as well: ogrinfo -ro -sql ...


2

Your 'ENGL_NAME' shouldn't be abbreviated at all (less than 10 characters), but writeOGR has its own will, it seems. Instead of writeOGR(shp, "PolygonsV2", speciesname, driver="ESRI Shapefile") you might try currdir <- getwd() #store your current working directory setwd(paste(currdir,"PolygonsV2",sep="/")) #switch to your desired folder ...


2

Just figured out a way to do this in R, inspired by the link I posted in the comments (which uses outdated functionalities but the right packages); as an example, I'll use the "Southwesternmost" counties in Michigan (shapefile here--Allegan, Berrien, Cass, Kalamazoo, St. Joseph, and Van Buren counties) So, ordering counties alphabetically, the (length>0, ...


2

Ok I could get it done with numpy instead of looping. I've done this way: def useNumpy(pt1, pt2): #The input are two lists of points defining the two polygons of interest #Shorten some function names that will be repeated nr = numpy.roll nt = numpy.tile #Give lists as array arr1 = numpy.array(pt1) arr2 = numpy.array(pt2) A0 ...


2

You have many solutions but you choose to use ogr and the GeoJSON format (you could have chosen Fiona and shapely or one of the many modules to parse gpx files: Pypi:gpx or to parse XML files). 1) Is there a more direct or idiomatic way to create this rectangle? ogr is is verbose, look at Python GDAL/OGR Cookbook With the GeoJSON format: import json ...


2

Turns out it is as easy as this: MIMETYPE "application/json; subtype=geojson; charset=utf-8" Without the setting, no content encoding is returned by the server.


2

Stepping back, if you have data in PostGIS and you want to render it with Mapnik you really should be using Mapnik's native PostGIS support - aka mapnik.PostGIS datasource in python (https://github.com/mapnik/mapnik/wiki/PostGIS). That will be faster and more efficient than going through OGR. That said, if you have a good reason for connecting to PostGIS ...


2

Depending on your version of mapnik, yes it's possible. Mapnik 2.1 introduced the Python Plugin as a data source. You would need to subclass mapnik.PythonDatasource and implement a features method that based on an incoming mapnik.Query object will return mapnik.PythonDatasource.wkb_features objects. For example: import mapnik import ogr class ...


2

You on Linux? If so, switch all the single quotes for double quotes and vice versa like this: ogr2ogr -f GeoJson -where 'neighborhood IN ("Lower East Side", "Greenwich Village", "Columbia St", "Financial District", "Flatiron District", "Williamsburg", "West Village", "Central Park", "Upper West Side", "Navy Yard", "Gramercy", "Stuyvesant Town", "Upper East ...


2

This code works for creating an OGRDataSource. Special thanks to Michael Miles-Stimson for getting me on the right track. const char* path = "C:/Test/test.shp"; OGRRegisterAll(); OGRDataSource *hDS; OGRSFDriver *driver; OGRSFDriverRegistrar *registrar = OGRSFDriverRegistrar::GetRegistrar(); driver = registrar->GetDriverByName("ESRI Shapefile"); hDS = ...


2

do you mean a line geometry? something like: import numpy as np from osgeo import ogr v=np.array([[[[-10885205.10690245,3525662.26531131], [-10885369.01690829,3525374.0741439], [-10885424.08206484,3525269.35422813], [-10885439.88685361,3525278.66051977], [-10885600.07792469,3525092.02164548], ...


2

Don't expect OGR to be Pythonic. In fact, there is a large list of gotchas documented to describe these. Essentially, all of the variables are references to objects described in C++, so they often don't make sense with workflows in Python. Modules like copy only work on Python objects, which is why your example doesn't work as expected. You need to keep ...


2

In order to use this data in PostGIS or another system you need to convert the start-center-end arcs into start-midpoint-end arcs. This is actually pretty easy. In vector terms, you subtract the center from the start- and end-points. Now you can get the mid-point of the arc by adding the start- and end-points together, then normalizing to the radius of ...


2

That solved my problem. Further information and a tutorial can be found here http://gdal.org/1.11/ogr/ogr_apitut.html A C++ version for GDAL 1.11: #include <GDAL/ogrsf_frmts.h> int main() { OGRRegisterAll(); OGRDataSource *poDS; poDS = OGRSFDriverRegistrar::Open( "data.shp", FALSE); }


2

You can test it with gdalsrsinfo http://www.gdal.org/gdalsrsinfo.html gdalsrsinfo epsg:3068 PROJ.4 : '+proj=cass +lat_0=52.41864827777778 +lon_0=13.62720366666667 +x_0=40000 +y_0=10000 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.45 5,6.7 +units=m +no_defs ' OGC WKT : PROJCS["DHDN / Soldner Berlin", GEOGCS["DHDN", ...


2

A MultiPoint, by definition, comprises several points, so geom = '{"type": "MultiPoint", "coordinates": [[100.0, 0.0], [101.0, 1.0]]}' geom = ogr.CreateGeometryFromJson(geom) # iterate through the points for i in range(0, geom.GetGeometryCount()): g = geom.GetGeometryRef(i) print "%d,%d" % (g.GetX(), g.GetY()) 100,0 101,1 as explained in Iterate ...


2

If you only want the point coordinates, you don't need Shapely, simply use the appropriate key of the dictionary: for point in filter_list: print point['geometry']['coordinates'] (270977.604378, 153144.810665),... If you want a Shapely geometry, use the shape function of Shapely: from shapely.geometry import Point, shape for point in filter_list: ...


2

Your capital_pt is the coords attribute of the original capital shapely geometry object. In itself this is not a shapely geometry, rather a sequence of tuples of flots which are the point objects. Instead you should be using distance_between_pts = capital.distance(city_items) Given your case where you're trying to calculate distances for all points in a ...



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