New answers tagged

1

I started it all again: restored a previous database from backup; create extension postgis; create extension postgis_topology; ogr2ogr from console creates additional tables (like waypoints); after this I need to grant access in those additional tables to php/apache user.


0

Below is the answer to my problem. hope it helps someone def findCountryByBoundingBoxOffline(coords): log = logging.getLogger(__name__) bb = getBoundingBox(coords) path = 'borders' file = os.path.join(os.path.dirname(__file__),path, 'World_Country_Borders_KML.kml') #SpatialRef = osr.SpatialReference() ...


3

What is strange if the intersection result is a point ? The intersect predicate is Returns True if the boundary and interior of the object intersect in any way with those of the other. With a common point between the geometries, the intersects predicate returns TRUE because the boundary of the first geometry intersects the boundary of the second ...


1

Not a proper answer but does not fit into the comment box. However, at least PostGIS does return a point for two polygons which touch at one point. SELECT ST_AsText( ST_Intersection(ST_GeomFromText( 'POLYGON (( 140 360, 140 480, 220 480, 220 360, 140 360 ))'), ST_GeomFromText( 'POLYGON (( 220 260, 220 360, 300 360, 300 260, 220 260 ))'))); ...


1

Figuring out the GPX driver options for rgdal is headache-inducing. Writing a linestring as you've done here will cause it to write a route layer - if you write a multilinestring it should create a track layer. According to the documentation you should be able to make it be a track layer regardless using FORCE_GPX_TRACK=true but I've not been able to make ...


1

Actually it works, just need to wait few minutes after code exercise.


3

You don't need Gdal/Python here. It is easier to use Shapely 1) Transform the list of points to a shapely geometry list = [(0.0, 0.0), (0.0, 1.0), (0.0, 2.0), (0.0, 3.0), (0.0, 4.0), (0.0, 5.0), (0.0, 6.0), (0.0, 7.0), (0.0, 8.0), (0.0, 9.0), (1.0, 0.0), (1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0), (1.0, 5.0), (1.0, 6.0), (1.0, 7.0), (1.0, 8.0), (1.0, ...


2

You can do everything that is possible to do with SQL also with ogrinfo by utilizing the SQLite dialect http://www.gdal.org/ogr_sql_sqlite.html. Once you know that the problem turns into a question about how to make such a query and the answer can be found from Stackoverflow http://stackoverflow.com/questions/17194145/sql-count-based-on-column-value The ...


1

For your first question, I believe this previous asked question will help you: Why does this simple Python OGR code create an empty polygon? Here is a python GDAL cookbook (a few years old, but still applicable) GDAL/OGR python Cookbook: https://pcjericks.github.io/py-gdalogr-cookbook/ For your second question, I believe that you want to set a spatial ...


1

You can build a custom mercator projection centered approximately on the center of the swath. For example, use for swath 25: +proj=merc +lon_0=-140 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs In this projection, the swath is not broken by the dateline. You can create the polygon from the line. Then create a cut polygon ...


1

I would rewrite the swathe line generation process to start and finish in the same continuous longintudinal space. ie if a line started at 170° and finished at -170° I would rewrite the process to finish at 190° instead without wrapping at -180,180 Then you can make unbroken polygons between your lines. Then use a clip process to split the polygons at the ...


2

Be modern, use easiest modules for that (look at Problem intersecting shapefiles with OGR in Python) from shapely.geometry import shape import fiona with fiona.open("ReBu.shp") as input: meta = input.meta with fiona.open('result.shp', 'w',**meta) as output: for feature in input: if abs(shape(feature['geometry']).area - 39.8953) < ...


2

It's always risky to attempt equivalence operations on floating-point values. Even if the value displays "39.8953", it's likely to be stored as "39.8952999999999999997" or "39.895300000000000001". The only safe way to compare against floating-point values is perform range-checking: myQuery = 'Area_ha >= 39.89525 and Area_ha < 39.89535' Depending ...


0

I usually try to convert any files to Shapefiles via ogr2ogr and then load it to PostgreSQL, usingogr2ogrorshp2pgsql`. A little bit complex, but i save files for history and it's easy to see changes in attributes. ogr2ogr -f "ESRI Shapefile" your_shape.shp your_qpx.qpx -sql "SELECT somefieldfromgpx AS somefieldfromdatabase FROM your_qpx" note that ...



Top 50 recent answers are included