Hot answers tagged

5

Linux (ubuntu as you call it) uses a shell (probably bash in this case) which has a different syntax to windows so you want something like: for i in se70*.tif do gdalinfo $i done


5

You are trying to execute a command line utility from within Python. To do this you can use subprocess, which takes the command line arguments as a list of strings. import subprocess cmd = ["gdal_translate", "-of", "XYZ", "dataset.tif", "D:\Delft3D\Teste\GDAL\clc_raster.tif.xyz"] subprocess.call(cmd) edit: Additionally I'd like to point out that a ...


5

You have now a better way to do. Since RFC 59.1 : GDAL/OGR utilities as a library, you can use gdalwarp from Python directly without using any call to the command line utility but using really the function from Python. This solution is a bit "on the edge" as you need at the moment to use the latest GDAL version (version 2.1, in fact the master/trunk ...


5

I do not understand why beginners try to start with the GDAL/OGR Python bindings (not very "Pythonic" and difficult) when other easier alternatives are available. With your script, you need to know osgeo.ogr and the SQL dialect of SQLite. The solution proposed by Mike T is powerful but not "classic" and performs only the intersection of shapefiles. What ...


4

Rather than building as a plugin for an older release, you could download and use a GDAL 2.1 development snapshot build from GISinternals that contains the new Sentinel-2 driver. Both the MSVC 2012 and MSVC 2013 builds have it compiled in. gdalinfo --version GDAL 2.1.0dev, released 2015/99/99 gdalinfo --formats|findstr /I sentinel SAFE -raster- (rov): ...


4

Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where: # reclassification lista[np.where( lista < 200 )] = 1 lista[np.where((200 < lista) & (lista < 400)) ] = 2 lista[np.where((400 < lista) & (lista < 600)) ] = 3 lista[np.where((600 < lista) & (lista < 800)) ] = 4 ...


4

I don't understand your problem import osgeo.ogr as ogr import osgeo.osr as osr driver = ogr.GetDriverByName("ESRI Shapefile") srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # create the data source data_source = driver.CreateDataSource("test_ogr.shp") # create the layer layer = data_source.CreateLayer("test_ogr", srs, ogr.wkbPoint) # Add the field ...


3

The spatial reference of the tif is not defined. You will need to define it before proceeding with the clip. C:\Users\frogman\Downloads>gdalinfo Hillshade_Wrap.tif Driver: GTiff/GeoTIFF Files: Hillshade_Wrap.tif Size is 9600, 15120 **Coordinate System is ** Origin = (1773760.000000000000000,5897760.000000000000000) Pixel Size = ...


3

To convert an ASCII file with longitude, latitude and data value you may use a function like this: from osgeo import gdal def csv2tif(source, target): cvs = gdal.Open(source) if cvs is None: print 'ERROR: Unable to open %s' % source return geotiff = gdal.GetDriverByName("GTiff") if geotiff is None: print 'ERROR: ...


2

Sorry for the short answer but I suppose that simply adding an alpha channel might do the trick. At least I could make an output that is semi-transparent in QGIS with the following command. gdal_rasterize -of gtiff -ot byte -co alpha=yes -burn 255 -burn 0 -burn 0 -burn 100 -ts 500 500 -l test test.shp test.tif


2

If QGIS asks you for a CRS, it does not write that information into the file. You have to save the data into a new file, then QGIS writes the CRS information into it. Alternatively, use gdal_translate -a_srs to assign the CRS to the output file. If you need a different projection, use gdalwarp -s_srs -t_srs to do both things in one step.


2

you can use gdalbuildvrt instead of gdalmerge to create a virtual raster template (xml file) that will behave like a merged raster. then you can run gdal_translate with the -projwin option for your set of 256*256 tiles. For most of my applications, I would create a single output per face then use gdalbuildvrt for virtual tiles. Note that if you create one ...


2

You should be able to just use the sql argument in ogr2ogr. For instance, with the following polygon Shapefile with two features and two attributes: $>ogrinfo -so -al polygon.shp INFO: Open of `polygon.shp' using driver `ESRI Shapefile' successful. Layer name: polygon Geometry: Polygon Feature Count: 2 Extent: (-1.206294, -0.828671) - (0.727273, ...


2

You should exchange ymin and ymax in -te (target extent). -te 102.3375307079206 10.350077240783799 107.6323189079206 14.6874037407838


2

Your values aren't in 0,255 since they are UInt16. You can try rescaling to 0,255 (GDAL works it out by default from input min/max and output default 0,255): gdal_translate -b 1 -b 2 -b 3 -mask "none" "input.tif" "output.tif" -scale Note you can add params if the defaults aren't sensible: -scale [src_min src_max [dst_min dst_max]] ...


2

I'd check the following: first of all, I'd like to be sure that source data are OK, by comparing them with official/cadastral data I'd try to use the qgis openlayers plugin to check if the transformation to 3857/900913 works properly for both 23031 and 25831 Try to reproject using QGIS, just saving the shapefile with the target srid. If using QGIS works ...


2

Try clipping the polygons before using them (also, please try to provide complete code including library calls in the future): library(ggmap) library(rgdal) library(rgeos) library(ggplot2) URL <- ...


2

Whenever I need to process large amounts of data I run gdal from the command line. Check http://www.gdal.org for info, examples and a couple tutorials to get you comfortable using it from the command line.


2

Using the EPSG works: gdalwarp -t_srs EPSG:4269 RipBuf100_09.tif RipB uf100WGS_09.tif Creating output file that is 55104P x 24192L. Processing input file RipBuf100_09.tif. Using internal nodata values (e.g. -2.14748e+009) for image RipBuf100_09.tif. Copying nodata values from source RipBuf100_09.tif to destination ...


2

Solved using Programmatic raster-vector calculation I had some troubles using directly gdal.RasterizeLayer() on layers like propose here in the cookbook but it seems that using "MEM" data source and finally writing it on the disk is maybe better. My solution (which is one of many possible solutions) def rasterizer(shapePath, rasterPath, attribute, ...


2

You could use another option like: os.system("gdal_translate -of XYZ " + source_file + " " + out_file)


1

Put the proj4 string between quotes. Now the string gets split at parsing and GDAL thinks that the +ellps is a name of the source file.


1

shot in the dark, but maybe use the EPSG instead of the actual CRS label


1

The correct answer was given yesterday on IRC in the GDAL channel. The above *.nc file is a netCDF version 4 encapsulated into a HDF5 format. Under Linux the libraries support the version 4 of netCDF, so the file opens with this (specific) driver and the medatadata (as min/max values) are correctly displayed. Under Windows (at least with GDAL installed ...


1

I don't know why the order of drivers is different between Windows and Linux: Windows: gdalinfo --formats | findstr "netCDF HDF" HDF4 (ros): Hierarchical Data Format Release 4 HDF4Image (rw+): HDF4 Dataset HDF5 (ros): Hierarchical Data Format Release 5 HDF5Image (ro): HDF5 Dataset netCDF (rw+s): Network Common Data Format Linux: gdalinfo --formats|grep ...


1

So you just want to mosaic all the tiles for a given day? That's a perfect Job for GDALs VRTs. gdalbuildvrt mosaic_049.vrt 049*.hdf or from Python for all days import subprocess import glob for day in range(0, 365): day = str('%0.3d' % day) cmd = ["gdalbuildvrt", "mosaic_"+day+".vrt", glob.glob(day+"*.hdf"] subprocess.call(cmd) edit: ...


1

You are trying to adress a subdataset inside a HDF container directly. There are two ways you can do that with gdalinfo: Put the complete name of the subdataset in parantheses gdalinfo " HDF4_EOS:EOS_GRID:MYD13A2.A2015297.h16v05.005.2015314081208.hdf:MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI" Use the subdataset option gdalinfo -sd 1 ...


1

You have Georeferenced images that you want to Rectify. There doesn't seem to be a command line Tool that mimics the one on the Georeferencing Toolbar, Rectify. I found a script, by Rob, that reads the control Points from the .tif.aux.xml file and passes this info to the Warp command, outputing a warped image. Please see the script at the end of this post, ...


1

ArcGIS supports a number of methods to georeference an image. IMHO its a matter of taste, what methods are associated with georeferencing, and what with "undoing" projection or distortion. The simplest of them, the affine or 1st Order Polynomial supports translation, rotation, and scaling. This one is parameterized by the numbers saved in a "world file", ...


1

not the best solution but solved the problem for me. For anybody who wants to know how I did it: gdal_contour -i 1000 -off 0.05 input.grd output.shp -i = some very high number -off = minimum pixel value then: ogrinfo -al -so output.shp and parse the ogrinfo output to get the extent coords EDIT Better Solution (Python3 with Subprocess): gdal_calc.py ...



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