New answers tagged

0

Here is a auto color script based on AndreJ's code: ================================================================================ Usage: ------------- Custom color: python gdaldem.py input_tif.tif color.txt output_color.tif Auto color: python gdaldem.py input_tif.tif auto output_color.tif """ import subprocess import sys import os import tempfile import ...


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.


1

I think the following parameters requires a minimum value of 1: JPEGCOMPRESSION ZLEVEL PREDICTOR The others can be set to 0 so you could try running the following which works for me: processing.runalg('gdalogr:cliprasterbymasklayer', inputlayer, maskshape, "", False, False, False, 0, 0, 1, 1, 1, False, 0, False, "", "outputraster")


0

The actual remote path was like: private const String NetworkShapeFilePath = @"\\NetworkServer\FakePath\_vectordata.shp" and I had it written like: private const String NetworkShapeFilePath = @"\\NetworkServer\FakePath\vector_data.shp"


1

This answer isn't a single command, but I'll put it in to get the ball rolling. Use gdalwarp to resample the geoid grid, then gdal_calc.py to shift the original raster. gdalwarp -s_srs epsg:4326 -t_srs epsg:26910 -r cubic -tr 10 10 -tap HT2_0.gtx HT2_0_resampled.tif gdal_calc.py -A original.tif -B HT2_0_resampled.tif --calc="A+B" --outfile=shifted.tif If ...


1

I made a test with OSGeo4W version of GDAL gdalinfo --version GDAL 2.0.2, released 2016/01/26 I created a test point in OpenJUMP JML format with one decimal attribute. <?xml version='1.0' encoding='UTF-8'?> <JCSDataFile xmlns:gml="http://www.opengis.net/gml" xmlns:xsi="http://www.w3.org/2000/10/XMLSchema-instance" > ...


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


0

I believe it is possible but probably not very accurate. Multiple elevations do not have same RGB value if the ramp is continuous and not like having red at both extremes. I would probably start by reducing the number of colors and saving the result into a paletted image by using rgb2pct-py http://www.gdal.org/rgb2pct.html. Then I would vectorize the ...


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 ))'))); ...


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


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


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


5

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


1

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


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


0

I just solved this by placing the right files extracted from RPM. Here is where I get my 64bit RPM http://www.filewatcher.com/d/Mandriva/2010.1/x86_64/Sciences/Geosciences/gdal-python-1.7.2-1mdv2010.1.x86_64.rpm.220418.html Using the rpm extraction command provided by http://www.cyberciti.biz/tips/how-to-extract-an-rpm-package-without-installing-it.html I ...


1

I eventually gave up trying to load the large files - it does not seem to be possible with any tool. The other side of the coin was that the features I did manage to load took over 1 hr to do a simple STIntersects query. Solution: I cut up the large shapefiles using ogr2ogr based on county geometries and then uploaded them and associated them with the ...


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


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

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

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.


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!


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.


0

I had such problems with not only python gdal binds but with standard gdal apps: gdalwarp, gdal_translate and so on.. So the problem was solved after full reinstall of all gdal and python reminds in system with refreshing path and gdalpath enviroment vars. The fastest way: Install Anaconda https://www.continuum.io/downloads, Install GDAL core from ...


0

I was able to use the '-nlt linestring' to export a line feature class successfully. The error appears to have been embedded in the sql statement. I was joining two tables and a view. The view was originally written using the 'myview.*'. I changed the SQL so that added each variable individually. This seemed to work. table1.var1, table2.var2, ...


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


-1

I found solution: self.sxf_ds = gdal.OpenEx(self.sxf_path, gdal.OF_VERBOSE_ERROR, ['SXF']) As i understand: GDAL is C/C++ library and main docs are written for these languages. All methods(functions) and constants are linked via SWIG. gdal.py imports gdalconstants.py and in the second one on line 341 (gdal 2.0.1) there is declaration of "my" constant ...


0

As written here: "libpoppler itself must have been configured with --enable-xpdf-headers". Are you sure libpoppler has been compiled against the xpdf C++ headers?


0

This is an old question, but I've just run into what I think is a similar problem dividing two integer grids to get a percentage. The key for me was to convert them to float grids in the calculation. Just setting the output type didn't help. I am running GDAL 1.11.2. Example code is: gdal_calc.py \ -A "_t2.tiff" \ -B "_t3.tiff" \ ...


1

I think you may have the ogr2ogr syntax wrong. The last line suggests that ogr2ogr didn't expect that type of argument. I'm using GDAL 2.1.0dev, and this works for me: ogr2ogr --config CARTODB_API_KEY XXXX -f CartoDB CartoDB:wri-02 myfile.shp I think you just need to remove the quotes around the "CartoDB:Development", provided that your account is really ...


0

If you build gdal yourself, you should also build grass yourself. Using grass binaries which are linked against a different version of gdal than QGIS is certainly a very bad idea.


-1

Thank you everyone for your help. I have finally solved this issue, and I would like to explain the procedure to do the same. Clear previous installations of gdal in your Ubuntu system. $sudo apt-get remove gdal-bin $sudo pip uninstall gdal $sudo apt-get remove python-gdal Download GDAL source from this link . You can browse and download one of the ...


0

I found an answer for exactly this problem in the [gdal-dev] Mailing list: To avoid the fourth band being interpreted as ALPHA channel, the corresponding option ALPHA=NO can be specified for the GeoTIFF driver. To preserve the RGB order, apparently there is the option PHOTOMETRIC=RGB - for more read the mailing list. In my case the ALPHA=NO option was ...


0

You have specified degrees as units, but not added any information what coordinates your corner points have. Apart from that, +units=d is not listed with proj -lu and might be ignored. You may try gdal_translate with -a_ullr, or add GCP points with -gcp.


0

It requires an object that implements the Progress interface. Either do your own (implement setInfo and whatever else is required) or use a shipped one like processing.gui.MessageBarProgress.MessageBarProgress. from processing.gui.MessageBarProgress import MessageBarProgress mp = MessageBarProgress() GdalUtils.runGdal(command, mp)


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


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

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


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

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


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)


0

gdal_fillnodata Fill raster regions by interpolation from edges. gdal_fillnodata.py [-q] [-md max_distance] [-si smooth_iterations] [-o name=value] [-b band] srcfile [-nomask] [-mask filename] [-of format] [dstfile]


0

Hey if you are using ArcGIS, below link will give you code developed by me which does the same thing that MRT tool does. Modify the paths of your data input and shape file for cliping the desired area. https://github.com/HGIS4YOU/MODIS_DATA_NDVI_Arcpy-code.git



Top 50 recent answers are included