Tag Info

New answers tagged

0

The title ('Extract raster values ...') is not reflecting your code instructions because you are mixing vector and raster layers. Assuming that the answer is only for raster layers, I think that the best way to read in them (for any purpose with Python/GDAL) is by using a scanline and the unpack struct function. The code is more compact, the control is more ...


0

Sorry for the code format. I posted again: import os, sys try: from osgeo import ogr, gdal from osgeo.gdalconst import * os.chdir('/home/digd/Desktop/puntos') except ImportError: import ogr, gdal from gdalconst import * os.chdir('/home/digd/Desktop/puntos') driver = ogr.GetDriverByName('ESRI Shapefile') shp = ...


0

Ogr2ogr easily handle your data but shapefile format does not. Ogrinfo reveals that your data contains one geometry which is of type Geometrycollection. Collection is made by combining one linestring and one point. GEOMETRYCOLLECTION ( LINESTRING ( 120.9954464 6.040449 0, 120.9953874 6.0402356 0, 120.9957951 6.0401396 0, ...


1

This can be found from gdalinfo, except it has three forms: If rotation / shear coefficients (adfGeoTransform[2] and adfGeoTransform[4]) are zero, the output is simplified: Origin = (%.15f,%.15f) % adfGeoTransform[0], adfGeoTransform[3] Pixel Size = (%.15f,%.15f) % adfGeoTransform[1], adfGeoTransform[5] If rotation / shear coefficients are non-zero, ...


0

In C# you can access GDAL provided you have the libraries installed. Reference gdal_csharp.dll. This tutorial helps a little. using OSGeo.GDAL; // later inside the class OSGeo.GDAL.Gdal.AllRegister(); // IMPORTANT or nothing will work // open the ECW dataset in GDAL OSGeo.GDAL.Dataset pInDS = Gdal.Open(ECW_Path, OSGeo.GDAL.Access.GA_ReadOnly); // get the ...


1

Using gdal_translate you have to set at least two arguments: src_dataset - yours SPOT NDVI HDF file dst_dataset - output TIFF filename, for example spot_ndvi_250m.tif (you're missing this one, that's why you've got an "No target dataset specified" error. Then in order to resample your image into 250m pixel using Nearest Neighbour you should add another ...


1

gdalwarp -geoloc allows you to use the complete 2d-array of latlon as georeference. With that, you can use any target CRS to reproject your data to a commonly used projection. See my answer to this question for an example: How to match a raster NetCDF data with a vector layer in QGIS?


2

The tool works well. I modified slightly the code for producing only 3 raster (10, 20 and 30 Pixel Size) with a point shapefile with 5 features. import processing RasterPixelSize = (10.0,20.0,30.0) i=1 for rps in RasterPixelSize: processing.runalg("gdalogr:rasterize",\ "/home/zeito/pyqgis_data/point.shp",\ "id",None,1,rps,rps,6, ...


0

-r bilinear Did your hillshade worked if reprojected using gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:3857 -r bilinear input.tif reproj.tif ? It didn't for me. Resizing After some tests : data => hillshade : fine (good) data => resizing => hillshade : stripped (bad) data => resizing => reprojected => hillshade : stripped (bad) Commands for ...


0

The answer of SaultDon is right. You have to install libhdf4-alt-dev. HDF4 contains an old version of netcdf which is disabled in this package. This is disabled in that package because it clashes with more recent versions of netCDF. Since GDAL uses the headers of the recent version of CDF you have to use the -alt- version. This will not disable netCDF ...


1

The scale option is only needed for "unprojected" rasters with units in degrees. It is the ratio of height units used in the DEM (typically meters) to distance units in degrees. So for Mercator or other projected rasters, you can ignore the scale option, or use the default of 1 (no scale).


2

The part that is going wrong is where you assumed that the coordinates in 3857 are degrees (like -5.8, 41). The units of 3857 are metres. So you're asking for something tiny, off the map (near the origin). Lets look at a conversion, using pyproj: from pyproj import Proj, transform inProj = Proj(init='epsg:4326') outProj = Proj(init='epsg:3857') x1,y1 = ...


0

I'm using GDAL Version 1.11 in C++, that said my answer may not be relevant depending on the version you are using. The the Create() method is not supported by the PNG driver (I'm not familiar with the JPEG driver), so you must first create a GDALDataset with a different driver the use a the CreateCopy method from a PNG driver to get your PNG image. Also, ...


0

You could do this easily enough in C++. The way I would do it is to send the path of the file you want to crop, the path of the output/cropped raster, the top left coordinates of the cropped raster and the width and height of the cropped raster. An outline for your code could look something like this. . . void crop(const char *inputPath, const char ...


0

I've done this way... it is a little hard-coded but if nothing changes with the gdalinfo, it will work for UTM projected images! imagefile= /pathto/image p= subprocess.Popen(["gdalinfo", "%s"%imagefile], stdout=subprocess.PIPE) out,err= p.communicate() ul= out[out.find("Upper Left")+15:out.find("Upper Left")+38] lr= out[out.find("Lower ...


0

I suppose that you are happy with your script but you may find use for the Spatialite SQL dialect in the future. Let's take a shapefile with one multipolygon to start with. ogrinfo multipolygon.shp -al INFO: Open of `multipolygon.shp' using driver `ESRI Shapefile' successful. Layer name: multipolygon Geometry: Polygon Feature Count: 1 Extent: ...


0

Here is the code I used to get this to work. This is VB.NET code: Public Function SchoolDistrictBoundary() As List(Of List(Of Coordinate)) Dim district As Feature = schoolDistrict() Dim result = New List(Of List(Of Coordinate)) If district IsNot Nothing Then Dim geometry = district.GetGeometryRef() Dim type = geometry.GetGeometryType() ...


1

It appears that you are missing the '-l dgm-10-epsg-32633' before your vrt-file. In essence, it doesn't know which layer of the vrt-file to use (even if there only is one layer). See the example in the man-page of gdal_grid: gdal.org/gdal_grid.html


4

Python does not contain a separate library for GDAL but rather bindings to access the GDAL libraries. As GDAL is updated so are the affected bindings, ensuring that you'll have full access to GDALs functionality from Python. If you upgrade GDAL you will also upgrade the Python bindings (if the version you upgrade to supports Python bindings). If you build ...


0

If you don't specify the target resolution in gdalwarp, it makes an educated guess of what could make sense. In most cases, it uses an average of the x and y resolution of the source data. But if you reproject from meters to degrees, this might go wrong. So you better set the target resolution or size with tror -ts (but not both!) explicitely each time. I ...


0

Reading this "Georeferencing from GeoTIFF is supported in the form of one tiepoint and pixel size, a transformation matrix, or a list of GCPs." I think is not possibile to write lon/lat for each image pixel like a layer...


1

It's telling you it doesn't like your values as int, they round to 0 and you can't have a 0 row, 0 column raster. That's because you're supplying the cell size and not the rows and columns. When you create a dataset it needs to know how many rows and columns it is, you then tell it how big the cells are in the geotransform: # 1. Define pixel_size and NoData ...


0

I would suggest using the MODIS reprojection tool (MRT - https://lpdaac.usgs.gov/tools/modis_reprojection_tool) to project and convert the data from HDF to TIFF. It's free (just need to create an account with NASA) and you'll be sure you're data is being transformed properly. Then you can work with your TIFF in R.


1

On Debian, you'll have to install libhdf4-alt-dev (not libhd4-dev) and then re-compile gdal. But there's a conflict with NETCDF if those drivers are enabled, so it's important to disable netcdf during the compile: ./configure --disable-netcdf --disable-fortran And it should find the hdf4 libs automatically. Reference: https://trac.osgeo.org/gdal/wiki/HDF ...


0

OK, it turns out there were two lines wrong in my code. The first is the line which selects the points iteratively from the pandas dataframe which should be: c = pts.ix[(pts.XX == XX) & (pts.YY == YY)] Before it was selecting point from the next grid cell and attributing them to the previous cell. The second error was in the geotransform line which ...


1

That solved my problem. Further information and a tutorial can be found here http://gdal.org/1.11/ogr/ogr_apitut.html A C++ version for GDAL 1.11: #include <GDAL/ogrsf_frmts.h> int main() { OGRRegisterAll(); OGRDataSource *poDS; poDS = OGRSFDriverRegistrar::Open( "data.shp", FALSE); }


1

If you look at this and this, they say that both GDALOpenEx() and GDAL_OF_VECTORS were (will be?) introduced in GDAL 2.0. GDAL 2.0 seems to be still under development. In case that you are able to compile it, you can find the source code here. In older versions you would use OGROpen to read a Shapefile.


0

If you want another file extension, you have to tweak the variable like this way: for %%N in (D:\Karten\gdal\gdal2tiles\NL25\*.tif) DO gdal_translate -of vrt -expand rgba %%N D:\Karten\gdal\gdal2tiles\NL25\%%~nN.vrt See http://ss64.com/nt/syntax-args.html for details on that.


2

You're expecting a 8-bit file, but you aren't exactly telling gdal that you want that. You need to add '-ot Byte'. See this help page for more info. Your command becomes: for %%f in (*.img) do gdal_translate -ot Byte -of ENVI -scale -1 1 1 255 -a_nodata 0 -co INTERLEAVE=BSQ "%%f" test/%%f.dat


0

You can always inspect the GDAL utilities for how they do it. For example, if you just want to "chip out" a section of an image, you could use gdal_translate -srcwin. The simplest way would just be to invoke gdal_translate (e.g. with a system() call, or maybe a CreateProcess() call if you're on Windows). Otherwise, you'd could read the code for ...


0

Assuming that GDAL was installed in your Windows system with a procedure similar to this: http://cartometric.com/blog/2011/10/17/install-gdal-on-windows/ precisely one way to confirm successfully installation of GDAL is appealing to gdalinfo --version in a MS-DOS Console. By the way, I am using GDAL 1.10.1. When I used a similar version of your command, ...


1

It works for me with a two-step-transformation: gdalwarp -overwrite -t_srs "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-100 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -of GTiff D:\Download\NEU\ster-wgs84.tif D:/Download/NEU/ster-100-transformed.tif gdalwarp -overwrite -wo SAMPLE_STEPS=1000 -t_srs "+proj=stere +lat_0=90 +lat_ts=60 ...


1

The ASCII raster needs a projection data file just as an ESRI shapefile. You can test it with one GTiff raster in QGIS using the Raster -> Conversion -> Translate (Convert Format) menu option where Arc/Info ASCII Grid format is chosen for output raster. You will see the *.prj associate file in the raster directory when conversion is over. On the other ...


1

Try with larger number of sample steps -wo SAMPLE_STEPS=1000 Reprojecting Natural Earth raster into natearth projection with gdalwarp causing black squares? Or with lower error threshold -et 0.01 gdalwarp leaves horizontal artifacts regridding from EASE-Grid (laea) to Polarstero (stere)


1

I've caught my error. May be it will be helpful for somebody. CreateCopy() makes read-only Dataset for PNG and JPG, so RasterIO() doesn't write anything into dataset's bands. We have to call RasterIO for the memory Dataset. I've made tests for PNG, JPEG and TIFF, it works. Create() call for "MEM" dataset doesn't allocate memory for the raster. Memory ...


1

One strategy is to use a nan_to_num function from Numpy, however it always uses 0. for nan, and there are unfortunately no parameters to replace nans with custom values. $ gdal_calc.py -A nan.tif --outfile=result.tif --calc="nan_to_num(A)"


1

This might be totally tedious but if you convert to Esri ASC (GDAL_Translate -of AAIGRID) which is a text based format and then open in a really good text editor like Notepad++ or VI/VIm, or write a python script, you can replace bad values using find & replace. Then save & close the file and translate to a GeoTiff


1

You can use the flag --NoDataValue=NODATAVALUE to replace NoDataValue. See the third example here


1

Like mentioned before, you might want to try OSSIM. It is a series of commandline programs and one or two GUI's. I am currently trying that out. There is no option to use a cutline as far as I know. Information about OSSIM: http://trac.osgeo.org/ossim Installers can be found at: http://download.osgeo.org/ossim/installers The current GUI is called ...


1

Not all band types are supported with all format. If you look at the format support pages for jpg and png, 8 or 16 bit are needed. Float values could be simply rounded so that it does not results in an error, but gives some zeros or NoData. You should try to use "rescale" to fit in a Byte. JPEG files are created using the "JPEG" driver code. Only Byte ...


0

As user30184 said, in Python, the process would be creating a memory raster of the same dimensions (layers and layer extension), and executing the CreateCopy after that: driver = gdal.GetDriverByName( 'MEM' ) driver2 = gdal.GetDriverByName( 'PNG' ) ds = driver.Create( '', 255, 255, 1, gdal.GDT_Int32) ds2 = driver.CreateCopy('/tmp/out.png', ds, 0) Then, ...


0

You should use a scanline and the unpack struct function with Python/GDAL. The scanline/struct method depends on fmttypes and their values can be supplied in a dictionary. For testing my answer I prepared a little raster of 29 rows x 29 columns and selected a region of 9 x 9 pixels in the arbitrary position (xoff,yoff) = (5,7). I registered manually the ...


0

I think the best way to read in a raster for any purpose with Python/GDAL is by using a scanline and the unpack struct function. The code is more compact, the control is more effective and the execution time is faster than the one with 'ReadAsArray'. The scanline/struct method depends on fmttypes and their values can be supplied in a dictionary. In the next ...


3

I guess you have gdal and the bindings installed, and some coding ability, so I'll just provide an outline: import org.gdal.gdal.gdal; import org.gdal.gdal.Band; import org.gdal.gdal.Dataset; ... Dataset dataset = gdal.Open(filename); Band band = dataset.GetRasterBand(1); ... // Do some band operation, like band.ReadRaster() to get the data, whatever you ...


3

Regarding your Post Scriptum: You already use the right projection. CRS:84 is synonymous to the EPSG-Code: 4326, which is the Coordinate-Reference-System known as: WSG84. You can always use EPSG-Codes instead of CRS-Names. EPSG-Codes are not so ambiguous in my opinion. In your Example, the transform-srs-parameter (-t_srs) would look like this: -t_srs ...


3

The clue is in the name :-) You downloaded boundary lines when infact you wanted to use States and Provinces which is a polygon layer.


1

If you want to write a lasfile with SRS information, you must create your liblas::Writer using a liblas::Header that has a defined SRS. Change your code to this: ofstream ofs; liblas::Header header; liblas::SpatialReference srs; srs.SetFromUserInput("EPSG:4326"); header.SetSRS(srs); header.SetPointRecordsCount(1); ofs.open("test.las", ios::out | ...


0

After a day of struggling with GDAL, I can confirm the methodology works and runs quickly. I was working with ice chart E00 files from CIS which I had converted to arcinfo coverage using avcimport. I have posted my prototype code below to help others. The E00 was converted to an arcinfo coverage called "test". The script must be run from the same directory ...


1

If you want to preserve the size of your raster, you have to specify that in the gdalwarp command line. Otherwise gdalwarp tries an own guess, based on the average of x and y source resolution preserving the extent, leading to square raster cells. It is pure coincidence that the pixel and line values are swapped in your case. You can force to preserve the ...


2

The OSGEO4W setup provides you with GDAL executables in Advanced install, including python bindings for python 2.7 and python 2.7 to run them. Gisinternals is another source, but not python.



Top 50 recent answers are included