Hot answers tagged

5

I usually look at http://epsg.io or http://spatialreference.org for all my WKT needs.


4

There is no command GDAL. GDAL is a library with many different tools within it. For a full list see: http://www.gdal.org/gdal_utilities.html If you just want the raster info, try: gdalinfo NDVI.tif


4

There is a configuration switch for this in GDAL 2.1 onwards: http://www.gdal.org/drv_geojson.html It defaults to 16 decimal places. So upgrading to GDAL 2.1 should solve your issue. If not you can try a higher number of decimal for floating point numbers: -lco SIGNIFICANT_FIGURES=17


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


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


3

GDAL has methods for that. From http://www.gdal.org/classOGRSpatialReference.html: "contains methods for converting between this object organization and well known text (WKT) format". ExportToWkt is probably what you need "Convert this SRS into WKT format. Note that the returned WKT string should be freed with OGRFree() or CPLFree() when no longer ...


3

Just my $0.02, there may be better ways. arcpy.env.overwriteOutput() Handle this yourself, ogr.Driver.DeleteDataSource() and ogr.DataSource.DeleteLayer() can handle this. You can use OGR_TRUNCATE, but this appears to be at the layer level arcpy.MakeFeatureLayer_management() If you have writable ogr.DataSource, then ogr.DataSource.CreateLayer(...) ...


2

I can't help with C# but I must warn you: what you believe to work correctly from the command line does not because you have comma as decimal separator. See the following examples: gdaltransform -s_srs epsg:4326 -t_srs epsg:3857 11,4110511566842 -16,3912563064512 1224514.39872601 -1804722.76625729 0 11.4110511566842 -16.3912563064512 1270272.40417808 ...


2

Your grid currently lacks projection information. Hence the line Coordinate System is '' in the output. With proper projection information, the output would start (for example) like: Coordinate System is: PROJCS["UTM Zone 4, Northern Hemisphere", ... It is important to understand the difference between assigning a projection and reprojecting. Using the ...


2

I think that the direct equivalencies are impossible unless you have a person who knows ArcPy and GDAL. Moreover the GDAL/OGR Python bindings (osgeo (GDAL/OGR) are not very "Pythonic" and difficult. It exists other easier alternatives ( Fiona, Pyshp (shapefile), GeoPandas,shapely, rasterio, ...). The last module is compared to ArcPy in Comparing Map ...


2

You have to point your system to where the GEOS Framework is located. See this post: DYLD_LIBRARY_PATH="/Library/Frameworks/GEOS.framework/Versions/3/unix/lib" export DYLD_LIBRARY_PATH If it works, add it to your $PATH!


2

--- gdal --- use gdal_contour function. documentation here ---- ArcGIS version----- posted before the software was specified Assuming you are using ArcGIS: 1. Convert the GeoTiff into ESRI GRID (esri raster format) using raster to other formats function, or simply import it in arcGIS and use save (export data) and save it as ESRI GRID. 2. use the Contour ...


2

The next code only select the number 3 band (blue band) in a RGB raster and write it as blueband.tif. from osgeo import gdal, osr import os, struct import numpy as np layer = iface.activeLayer() provider = layer.dataProvider() path = provider.dataSourceUri() fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', ...


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

the test yields either 1 (true) or 0 (false). So you can make an "and" by multiplying the results, and the othe cases using 1-(condition). Here ismy suggestion B*(A==212)*(B==20) + A*(1-(A==212)*(B==20))


1

Have you tried: projected coordinate systems: https://developers.arcgis.com/javascript/jshelp/pcs.html geographic coordinate systems: https://developers.arcgis.com/javascript/jshelp/gcs.html


1

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


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

The Shapely distributions I am making for OS X (https://pypi.python.org/pypi/Shapely#downloads) have GEOS included and you won't have to think about library paths at all. If you're using Python 2.7, 3.4, or 3.5 and OS X 10.6+, pip install shapely is the best way to get it.


1

Converting one State Plane to the adjacent system is always on the safe side, unless you need sub-meter accuracy. Alternatively, you can use EPSG:102004 USA_Contiguous_Lambert_Conformal_Conic +proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs which is valid for the most of the United States.


1

If you download the GDAL source, you'll have the source for gdal_translate.cpp. Just peruse the code and duplicate it in C#. The translation is straight-forward.


1

http://www.bev.gv.at/portal/page?_pageid=713,1572954&_dad=portal&_schema=PORTAL notes that the Austrian DEM is in MGI Lambert projection. https://www.data.gv.at/katalog/dataset/b5de6975-417b-4320-afdb-eb2a9e2a1dbf provdies the data in a 3 GB large Geotif file. This should be correctly georeferenced to EPSG:31287, and is free of charge. Some ...


1

Before installing GDAL python binding, you must install GDAL on your system. apt-get install -y libgdal-dev After that, set variables and install binding CPLUS_INCLUDE_PATH=/usr/include/gdal \ C_INCLUDE_PATH=/usr/include/gdal \ pip install GDAL Version of the system GDAL should be higher than the version bindings.


1

You must understand the logic of the Python version of GDAL. With a raster with the same projection (EPSG:3857) from osgeo import gdal data_set = gdal.Open("test.tif") # Getting the Spatial Reference System (returned in Well Known Text(WKT) format) projInfo = data_set .GetProjection() projInfo 'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS ...


1

If the output from gdalinfo file_name.tif is to be believed, it looks like the C++ implementation assumed your datum was WGS84 and filled out the rest of the WKT for you, although there shouldn't be an entry for both PROJCS and GEOGCS. It's hard to tell for certain without digging into the gdal source code and libraries.


1

If you need to know what these tags mean, there are a couple of documents on remotesensing.org here for 34735 and here for 33550. These seem to cover the meaning of the fields. As for how to set them programatically, that depends on which language you're using. The Geotiff wiki on osgeo.org has links to the spec and libraries for the Geotiff format.


1

Update Sometimes the solution is simpler than you think. First, I proceed as described by Joseph. After carrying out the converter function I get a .geojson file. This can be further processed as follows: json.load(open(pathGeoJSON)) From the generated JSON file the coordinates can be extracted and a QgsFeature generated.


1

You could try using something like the following which was tested in the Python console where you can call the output of the clip algorithm as an input to the convertformat algorithm: clip_layer= processing.runalg("qgis:clip", layer2, layer1, None) processing.runandload('gdalogr:convertformat', clip_layer["OUTPUT"], 1, "", None)


1

You can't (without modifying the code) as the exceptions module was removed in python 3. Note: the package metadata specifies Python 2 only: Programming Language :: Python :: 2 :: Only I suggest you raise an issue requesting a port to Python 3. You could try downloading, modifying and building the code yourself, then installing without root to your user ...


1

To create a new raster with gdal_rasterize, either the -tr or -ts option must be used, and the -b option cannot be used. E.g. gdal_rasterize -burn -9999 -tr 0.001 0.001 PG:"..." -l auckland auckland.tif



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