Hot answers tagged

126

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


58

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

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


48

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


45

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


42

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.


38

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


35

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


34

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(ds): """ Return list of corner coordinates from a gdal Dataset """ ...


33

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


32

With PyQGIS ras = QgsRasterLayer("raster.tif") pixelSizeX= ras.rasterUnitsPerPixelX() pixelSizeY = ras.rasterUnitsPerPixelY() print pixelSizeX 2.11668210081 print pixelSizeY 2.11685012701 With GDAL from osgeo import gdal raster = gdal.Open('raster.tif') gt =raster.GetGeoTransform() print gt (258012.37107330866, 2.11668210080698, 0.0, 163176.6385398821, 0....


30

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


29

Using rasterio you could do import rasterio file_list = ['file1.tif', 'file2.tif', 'file3.tif'] # Read metadata of first file with rasterio.open(file_list[0]) as src0: meta = src0.meta # Update meta to reflect the number of layers meta.update(count = len(file_list)) # Read each layer and write it to stack with rasterio.open('stack.tif', 'w', **meta) ...


28

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


28

The Gdal_translate utility can be used. The documentation mentions: ...to convert raster data between different formats, potentially performing some operations like subsettings, resampling, and rescaling pixels in the process. It also has an option for bands, where you selects which bands you want to operate on. So if you want to export just the ...


28

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


28

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


27

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.


27

If you don' want the values above 255 to be cut, you need to scale them down. For that purpose gdal_translate provides the option -scale: From the Manual: -scale [src_min src_max [dst_min dst_max]]: Rescale the input pixels values from the range src_min to src_max to the range dst_min to dst_max. If omitted the output range is 0 to 255. If ...


26

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


25

You could create a virtual mosaic from all Tiff files: gdalbuildvrt mosaic.vrt c:\data\....\*.tif and convert it afterwards: gdal_translate -of GTiff -co "COMPRESS=JPEG" -co "PHOTOMETRIC=YCBCR" -co "TILED=YES" mosaic.vrt mosaic.tif Keep an eye on all the GDAL creation parameters to compress your mosaic and use gdaladdo to add overviews. More info here: ...


24

Here's some python code that does what you want, reading GDAL files that represent data at specific times and writing to a single NetCDF file that is CF-Compliant #!/usr/bin/env python ''' Convert a bunch of GDAL readable grids to a NetCDF Time Series. Here we read a bunch of files that have names like: /usgs/data0/prism/1890-1899/us_tmin_1895.01 /usgs/...


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 ogr2ogr utility supports a limited sql syntax. You can join your CSV to the shapefile using something like the following: ogr2ogr -sql "select inshape.*, joincsv.* from inshape left join 'joincsv.csv'.joincsv on inshape.GISJOIN = joincsv.GISJOIN" shape_join.shp inshape.shp


24

You can download GDAL 2.1 for Windows from GIS Internals. There is an installer and a portable version that doesn't require installation. GDAL 2.1 is available for Ubuntu 16.04 from the UbuntuGIS-Unstable PPA sudo add-apt-repository -y ppa:ubuntugis/ubuntugis-unstable sudo apt update sudo apt upgrade # if you already have gdal 1.11 installed sudo apt ...


22

Condensed procedure outlined in http://cartometric.com/blog/2011/10/17/install-gdal-on-windows/ for Windows 7, 32 Bits, to install GDAL PYTHON: 1) Install Python. I installed Python 2.7.9 from https://www.python.org/ 2) Install the GDAL binaries published by Tamas Szekeres. First, I launched IDLE (Python GUI) noting the following values: "MSC v.1500" and ...


22

ogrinfo -so -al yourshapefile.shp This will give you geometry type, number of features/shapes, bounding box corners, projection information, and the name of each attribute file as well as the datatype stored in those attributes.


22

The problem is that I was not creating a field to store the raster band. After digging through the gdal_polygonize.py file, I realized this is not automatically done when calling gdal.Polygonize, which instead uses the function found here. Here is the extra step needed to create a field and write a band to the field: newField = ogr.FieldDefn('MYFLD', ogr....


21

OGR has its own idiom for stdin, /vsistdin/. Use that as ogr2ogr's first argument (the dst_datasource_name) and you can pipe curl's output to it: curl "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_50m_admin_0_countries.geojson" | ogr2ogr -f "KML" countries.kml /vsistdin/


21

You can use rasterio to extract the raster values within a polygon as in GIS SE: GDAL python cut geotiff image with geojson file I use here a one band raster file and GeoPandas for the shapefile ( instead of Fiona) import rasterio from rasterio.mask import mask import geopandas as gpd shapefile = gpd.read_file("extraction.shp") # extract the geometry in ...


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