Hot answers tagged

112

selimnairb's answer is close but you wont have the headers unless you've installed libgdal-dev: sudo apt-get install libgdal-dev with that done, export CPLUS_INCLUDE_PATH=/usr/include/gdal export C_INCLUDE_PATH=/usr/include/gdal and then pip install GDAL The compilation ran to completion and I have GDAL in my virtual env. Phew! (edit 2018) Note: To ...


90

The recommended way to do this is using the pgsql2shp utility, which should be installed with PostGIS. Note that you must include the geometry column in the query. $ pgsql2shp -f <path to output shapefile> -h <hostname> -u <username> -P <password> databasename "<query>" Example (creates qds_cnt.shp in current directory): $ ...


64

You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module. An example: from osgeo import gdal,ogr import struct src_filename = '/tmp/test.tif' shp_filename ...


50

This can be done in far fewer lines of code src = gdal.Open(path goes here) ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform() lrx = ulx + (src.RasterXSize * xres) lry = uly + (src.RasterYSize * yres) ulx, uly is the upper left corner, lrx, lry is the lower right corner The osr library (part of gdal) can be used to transform the points to any ...


49

In QGIS 2.8.1 there's a tool in the Vector menu that can split a dataset based on an attribute. Look in Vector > Data Management Tools > Split Vector Layer... It's a basic tool that should work if you don't want to resort to plugins or tools in Processing; unless they offer additional functionality you need.


47

See the OGR Projections tutorial and the OGRSpatialReference class. In particular, the GetAttrValue method. Here's a worked example. from osgeo import gdal,osr ds=gdal.Open(r'SOMERASTER.TIF') prj=ds.GetProjection() print prj srs=osr.SpatialReference(wkt=prj) if srs.IsProjected: print srs.GetAttrValue('projcs') print srs.GetAttrValue('geogcs') For my ...


46

if you have python-gdal bindings: import numpy as np from osgeo import gdal ds = gdal.Open("mypic.tif") myarray = np.array(ds.GetRasterBand(1).ReadAsArray()) And you're done: myarray.shape (2610,4583) myarray.size 11961630 myarray array([[ nan, nan, nan, ..., 0.38068664, 0.37952521, 0.14506227], [ nan, nan, ...


44

OSGEO4W and all standalone QGIS installers come with a OSGEO4W Shell. Start that, and type gdalinfo --version and read the result. You may have different versions on the disk: Standalone, OSGEO4W and also from gisinternals if you want the latest GDAL build, but every package sets its environment so that it is using the version it was delivered with. ...


42

Nutshell Each set of 3 images below should be read such as "grey (band) + opacity (band) = transparent result". You can test these processes within minutes via the associated github hosted makefile. Process #3 is the one which I recommend, with a threshold between 170 (keeps strong shadows) and 220 (keeps all shadows). Process 3 provides the strongest ...


38

You need to use the -dstalpha option to gdalwarp e.g.: gdalwarp -cutline INPUT.shp -crop_to_cutline -dstalpha INPUT.tif OUTPUT.tif This will add an alpha band to the output tiff which masks out the area falling outside the cutline. A late answer, but hopefully it will help someone else with the same problem.


36

Τhe header files cannot be found for some reason. Maybe you you are operating inside a Virtual Enviroment or they are not where they should be for some reason. In any case you can specify the include dirs when installing gdal via pip. first download python's gdal : pip install --no-install GDAL in later versions of pip (>= 9.0.0) pip install --no-install ...


35

Grab a copy of the ogr2ogr Python port, which is distributed with the GDAL source code download or can be found here: http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr2ogr.py Once you import that into your code, you can use it like this: import ogr2ogr def main(): #note: main is expecting sys.argv, where the first argument is the script name ...


33

Since GDAL 2.1 (more info here), GDAL and OGR utilities can be used as library functions. For instance: from osgeo import gdal ds = gdal.Open('input.tif') ds = gdal.Translate('output.tif', ds, projWin = [-75.3, 5.5, -73.5, 3.7]) ds = None


32

Below is an example that I wrote for a workshop that utilizes the numpy and gdal Python modules. It reads data from one .tif file into a numpy array, does a reclass of the values in the array and then writes it back out to a .tif. From your explanation, it sounds like you might have succeeded in writing out a valid file, but you just need to symbolize it ...


32

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


32

It's a well known and longstanding issue that gdalwarp doesn't deal with compression well. The solution is to gdalwarp without compression then gdal_translate with compression. To avoid two lengthy processes, gdalwarp to VRT first, it's really quick, then gdal_translate with the -co compress=lzw option. i.e. $ gdalwarp -tap -tr 30 30 -t_srs "etc..." -of ...


31

The shapefile should have a .prj file which defines the projection. You can use it together with one of the following 3 options to get either the proj4 string, WKT definition or EPSG code. To get proj4 definition: If you have gdal installed on your system, you can use the gdalsrsinfo command line application to get the proj4 definition as the OGC WKT ...


30

Here's another way to do it without calling an external program. What this does is get the coordinates of the four corners from the geotransform and reproject them to lon/lat using osr.CoordinateTransformation. from osgeo import gdal,ogr,osr def GetExtent(gt,cols,rows): ''' Return list of corner coordinates from a geotransform @type gt: C{...


29

In pure Python, without using the subprocess module (os.system is deprecated) to call ogr2ogr or shp2pgsql, for example): 1) with ogr Append a shapefile to a postgis table using the GDAL/OGR Python interface 2) with ogr and psycopg2 from the book Python Geospatial Development (Eric Westra), Chapter 7, p.219 import os.path import psycopg2 import ...


27

This is a known issue as documented on KyngChaos https://www.kyngchaos.com/blog/2018/20180322_qgis_3.0.0-3 A workaround is available as follows Go to Settings ... Options... System ... Environment Enable "Use Custom Variables " First select "Prepend", under variable enter "PATH", under value enter "/Library/Frameworks/GDAL.framework/Programs:/...


26

See the GDAL API Tutorial. #Import gdal from osgeo import gdal #Open existing dataset src_ds = gdal.Open( src_filename ) #Open output format driver, see gdal_translate --formats for list format = "GTiff" driver = gdal.GetDriverByName( format ) #Output to new format dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 ) #Properly close the datasets to ...


26

1) individual shapefile: as in the comment, a shapefile has only one layer. If you want only the names of the fields from osgeo import ogr source = ogr.Open("a_shapefile.shp") layer = source.GetLayer() schema = [] ldefn = layer.GetLayerDefn() for n in range(ldefn.GetFieldCount()): fdefn = ldefn.GetFieldDefn(n) schema.append(fdefn.name) print schema [...


25

Only adding this because I tried using the kyng chaos tools, but on my Mac OS X machine I was able to very, very easily install this with Anaconda conda install gdal Posting in case anyone finds this again - I realize the original post is 3 years old.


25

You can use rasterio to interface with NumPy arrays. To read a raster to an array: import rasterio with rasterio.open('/path/to/raster.tif', 'r') as ds: arr = ds.read() # read all raster values print(arr.shape) # this is a 3D numpy array, with dimensions [band, row, col] This will read everything into a 3D numpy array arr, with dimensions [band, ...


25

One possible solution to your problem: Convert it into a ASCII Raster, documention for which is here. This should be fairly easy to do with python. So with your example data above, you'd end up with the following in a .asc file: ncols 4 nrows 4 xllcorner 20 yllcorner 8.5 cellsize 0.5 nodata_value -9999 0.1 0.2 0.3 0.4 0.2 0.3 0.4 0.5 0.3 0.4 0.5 0.6 0.4 0....


24

follow this way: Find your spatial reference code from here Learn your tif file upper left coordinate and lower right coordinate. use this command for making it coordinated: gdal_translate -of GTiff -a_ullr ullon ullat lrlon lrlat -a_srs EPSG:4269 input.tif output.tif use this command for mercator:(epsg:3857) gdalwarp -of GTiff -t_srs EPSG:3857 input....


24

I've finally hit upon this solution, which I gained from this discussion (http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html). I like it because I can go straight from a numpy array to a tif raster file. I'd be very grateful for comments that could improve on the solution. I'll post it here in case anyone else searches for a ...


24

You could explicitly set the output coordinate range using the target extent option to gdalwarp (ie. "-te -180 -90 180 90") but you can also use the CENTER_LONG configuration option to force rewrapping around a new central longitude. Something like this: gdalwarp -t_srs WGS84 ~/0_360.tif 180.tif -wo SOURCE_EXTRA=1000 \ --config CENTER_LONG 0 ...


24

Quoting Frank Warmerdam, the maintainer: I pronounce it "goodle". I had originally thought to call it the "Geospatial Object Oriented Data Abstraction Library" (GOODAL) to make the right sound obvious, but I was too lazy to type GOODAL all the time, so I dropped the OO part. Some folks might say I dropped it from more than the name. :-) Most ...


24

The second argument in Open specifies if the data can be updated (written to). Try: dataSource = driver.Open(fn,1)


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