Hot answers tagged

5

This can be achieved with the help of GDAL's Virtual Raster Format. With this you can essentially skip the step of creating one giant DEM. The VRT will be handled by GDAL like a giant, merged DEM but is just a small XML file containing the file paths for each tile as well as some metadata. This can then be fed to gdalwarp together with a bounding box or a ...


4

GeoTransform is an array that contains six numbers: X origin Pixel width Angle (with vertical axis) Y origin Pixel height Angle (with horizontal axis) X and Y are coordinates of the top left corner of your raster image. Each pixel has size, width and height and it is a value in meters or degrees (depends on your CRS) - this value is constant, not ...


2

Have you tried using an equal blocksize. I deal with raster data which is of the order of 200k x 200k pixels and quite sparse. A lot of benchmarking has yielded 256x256 pixels blocks as most efficient for our processes. This is all to do with how many disk seeks are required to retrieve a block. If the block is too large then it is harder to write it to disk ...


2

Provide the source DEM no data value as well in the command line. Try the following command: gdalwarp -t_srs EPSG:4326 -srcnodata {source no data value} -dstnodata -90909 -of GTiff Source.tif DestWGS84.tif This should help in fixing the zero elevation issue you have.


2

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. P.S. duplicate question


2

I think all you need to do is reproject it using: gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -t_srs EPSG:3857 newImage.tif image.tif and then tile it: gdal2tiles.py newImage.tif If your file is very large it make take a while.


2

There is no fault in the data or python script. QGIS's default for opening a raster dataset looks for a minimum and maximum value to provide band rendering. However, this dataset contains a single value. Within QGIS's layer properties style menu by setting both the Min and Max values of the dataset to the single dataset value and selecting "No enhancement"...


2

You can do something similar using gdal_calc.py, e.g.: gdal_calc.py -A dtm.tif --calc='((A>=100)*(A<=200))*A+((A<100)*0)+((A>200)*0)' --outfile=dtm_reclass.tif --NoDataValue=-32767 This calc expression would: Assign a value of 0 to all pixel values less than 100 ((A<100)*0) Assign a value of 0 to all pixel values greater than 200 ((A>...


2

You are looking for GDALRasterBand::RasterIO. For efficiency, if you are indexing multiple points, you will want to read data in blocks then index into the resulting array


2

As far as I can see in the code, this is not possible with the HFA driver (with the GTiff driver, you wouldn't have this issue). On Linux, creation of the whole file is instantaneous as the file is created in a sparse way, but on Windows, it requires extra code to enable sparse mode from what I can see in https://msdn.microsoft.com/en-us/library/windows/...


2

I suggest to build a VRT around the full datasource: gdal_translate -of VRT http://user:pass@path_to_file/w1001001.adf local_path.vrt and work on with that. GDAL tries to write its .aux.xml file next to the source file, and will fail on that with an external source.


2

You are overwriting your data 'matrix'. First you reclassify the negative to 1, then you do other steps... but at the end yo are doing : matrix[np.where((matrix > 1.5)) ] = 6 Which is reclassifying your already changed values to 6. That is why you have only 1 (negative) and 6 the other values (2,3,4,5) -> 6 I would suggest to copy your matrix like ...


2

I've found when something isn't particularly well documented in GDAL, that looking through their tests can be useful. The /vsis3 test module has some simple examples, though it doesn't have any examples of actually reading chunks. I've cobbled together the code below based on the test module, but I'm unable to test as GDAL /vsis3 requires credentials ...


2

I didn't knew geopandas, it is very easy to add a field, a great library indeed. You just need to read your file add it and save. dataSrc = gpd.read_file('my_shp.shp') dataSrc['new_field'] = 1 dataSrc.to_file('newfile.shp') Selecting is also straightforward: dataSrc = gpd.read_file('my_shp.shp') dataSrc[dataSrc['id']<4].to_file('new_shape.shp')


2

GeoPDF is a trade mark of TerraGo and only programs which are licensed by TerraGo can create GeoPDF (tm) files. There are two systems to encode georeferencing in PDF and they are documented in http://www.gdal.org/frmt_pdf.html. Alternatives are ISO32000 or OGC Best Practice. GDAL can write the both variants if desired with a creation option parameter: ...


2

GDAL has no problem with GeoJSON http://geojson.org/geojson-spec.html. The following GeoJSON encodes one multilinestring: { "type": "FeatureCollection", "features": [ { "type": "Feature", "geometry": {"type":"MultiLineString","coordinates":[[[11.20558631,46.48251782],[11.2058444,46.48280049]],[[11.20578705,46.48252192],[11.20596731,46.48275133]]],"crs":{"...


1

Unfortunately, it turns out at this time, GDAL does NOT support OPENCL speed up on orthorectification. In fact, the only speedup I have been able to find is that compiling with OpenCL does speed up resampling slightly. This would be a good feature for GDAL to implement at a future date.


1

Try using an XML file to store the WMS info in, more details are at the GDAL WMS documentation. Here's an example WMS XML file to retrieve data from Mapzen's Elevation API: <GDAL_WMS> <Service name="TMS"> <ServerUrl>https://s3.amazonaws.com/elevation-tiles-prod/geotiff/${z}/${x}/${y}.tif</ServerUrl> </Service> <...


1

I would try to add a low resolution DEM (e.g. SRTM 30m) to fill the gaps where your high resolution DEM is not available. The problem with a crop before orthorectification is that the RPC link image coordinates with ground, so if your image is not georeferenced and you crop it, this relationship will be detroyed.


1

Personally, I've started using rasterio's windowed read/write However, you can use x and y offsets when writing your array to the output dataset. dataset.GetRasterBand(1).WriteArray(array, xoff, yoff)


1

Next code is similar to the cookbook code for testing the OGR 'Transform' method (from EPSG:32612 to EPSG:4326). It also includes 'AutoIdentifyEPSG' and 'GetAuthorityCode' methods for determining EPSG code for projection in my GeoJSON layer (EPSG:32612). from osgeo import ogr, osr import os driver = ogr.GetDriverByName('GeoJSON') # get the input layer ...


1

You can use the Zonal Statistics tool of QGIS (also available through the Processing Toolbox).. Giving the raster in input and choosing the vector with all your polygons, the algorithm will write for each feature (so for each polygon) some raster statistic, one of them is the max value.


1

When you set the Keep resoltuion option to True, you enable two parameters (according to the documentation for gdalwarp which is used in this function): -tr - Sets the output file resolution in target CRS georeferenced units. -tap - (GDAL >= 1.8.0) (target aligned pixels) align the coordinates of the extent of the output file to the values of -tr, such ...


1

I suppose that your image is some of the the 4-band products by Airbus DS: http://www.intelligence-airbusds.com/en/4951-which-spectral-mode-do-i-choose Gdal2tiles is made for splitting common, visual images into png tiles. Such images use 8 bits per band and they have one band (greyscale), 3 bands (red-green-blue) of 4 bands (reg-green-blue + alpha). I ...


1

I'm not sure if you can do this with any of the gdal cli tools, but I wrote something in python which accomplishes it: from osgeo import gdal from osgeo.gdalconst import GDT_Float32 import sys import numpy as np def fix_dem_nodata(raster_input, raster_output, nodata=0, threshold=-900): try: in_data, out_data = None, None # open ...


1

Create an ogr 'Memory' layer to rasterize from. Here's an example from the GDAL/OGR test suite: # Create a memory layer to rasterize from. rast_ogr_ds = \ ogr.GetDriverByName('Memory').CreateDataSource( 'wrk' ) rast_mem_lyr = rast_ogr_ds.CreateLayer( 'poly', srs=sr ) # Add a polygon. wkt_geom = 'POLYGON((1020 1030,1020 1045,1050 1045,1050 1030,...


1

The accepted answer is probably superior, but I just chose a quick and dirty method to merge files. First I build a index file using gdaltileindex. I then use gdal tools to returns a list of files that the bounding box intersects with. I then use gdal_merge.py to merge them the intersecting DEMs into one large dem and used the crop parameter of gdalmerge.py ...


1

1) I just don't know how to convert .geojson file to .shp. This is a one of the bases of ogr Python. If you have a geometry, it is very easy to convert it to a shapefile # geojson is GeoJson Polygon from osgeo import ogr output = "geojson.shp" driver = ogr.GetDriverByName('ESRI Shapefile') if os.path.exists(output): driver.DeleteDataSource(output) ...


1

WKT is only a text markup language for representing vector geometry, therefore you cannot convert 3D WKT to 2D (= text to text), you need to change the geometry The problem with ogr is that it seems that all the geometries are 3D by defaut 1) Creation of a 2D point with ogr point = ogr.Geometry(ogr.wkbPoint) # 2D point point.AddPoint(1198054.34, 648493.09)...


1

There are now Python modules easier to use for that, as rasterio Rasterio employs GDAL to read and writes files using GeoTIFF and many other formats. Its API uses familiar Python and SciPy interfaces and idioms like context managers, iterators, and ndarrays. Therefore from Masking raster with a polygon feature in Rasterio Cookbook import rasterio ...



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