New answers tagged

1

You can accomplish this using gdal and numpy. First you have to read all your rasters as numpy arrays and then stack them to create a single 3D array. Then you can mask this array to exclude values outside of the 0 - 30,000 range and then you can get the maximum value for each pixel. Finally, you have to fill the remaining masked values (pixels where the ...


1

Much simpler and intuitive if you use fiona and geopandas import fiona import geopandas as gpd # Get all the layers from the .gdb file layers = fiona.listlayers(gdb_file) for layer in layers: gdf = gpd.read_file(gdb_file,layer=layer) # Do stuff with the gdf Note: fiona uses gdal and geopandas uses fiona See also Reading the names of ...


2

SELECT ST_VALUE (rast,1,(ST_GeomFromText('POINT(319006 212230)',27700))) FROM public.country_rast WHERE ST_Intersects(rast, ST_GeomFromText('POINT(319006 212230)',27700))


2

What you did does edit the bounds but by the same you have destroyed the georeferencing. The original image is georeferenced with 4 ground control points and if you look at the coordinates you will notice that the image was not north-up but rotated/skewed. The right thing to do is to warp the image into a north-up image with gdalwarp gdalwarp -of GTiff -...


1

A raster is a set of cells that form a grid; each cell has a value. When you reproject a raster, you are re-drawing the grid to be aligned with a new projection. So, in the below figure, your original raster grid is shown in blue, and the reprojected grid is shown in red. Right away you can see a problem--the grids do not align. So, for example, in the ...


0

Found the answer. I used the -a_srs argument in gdal_transform and supplied minimum/maximum lat/longs. This Assigns/override the georeferenced bounds of the output file.


0

I tried an alternative approach where, rather than setting a spatial filter on my mask layer and rasterising that result, I decided to use ds.ExecuteSQL on my mask data source object. I constructed a spatialite query to find the intersection of my mask_ds geometries with the bounding box of my rast_ds. I then used this result layer to pass to the gdal....


2

When you iterate over the layer it essentially calls lyr.GetNextFeature(). After the last feature the pointer remains at the end which is why when you iterate over the layer the second time it doesn't find any features. If you plan on iterating over the layer more than once you need to make sure you call lyr.ResetReading() to reset the pointer to the first ...


1

The iteration over the features seem to consume these features like in an interator. Thus omitting the iteration for printing will solve the problem.


1

From the command line, execute this command inside the directory of GeoTIFF rasters you wish to stack (like a sandwich) together: gdal_merge.py -o your_depth_stack.tif -ps x y -separate *.tif This will stack all images inside the directory. If you only wish to stack a certain subset replace *.tif with the desired images: image1.tif image2.tif image3.tif ...


0

sometimes if you have used the ' ' in the filename like 'rainfall UK', this error would be happen. But you try to change it as 'rainfall_UK', it would be fine


1

Assuming your rasters are integer-type, you can create a ColorTable, specify the color for each value using the SetColorEntry() method and then apply the ColorTable to the raster using the SetRasterColorTable() method to the individual band. The SetColorEntry(pixel_val, (r, g, b)) method takes two arguments, where the first one is the pixel value and the ...


1

Reprojecting rasters of this type will always look "odd", if you need a projected GB base map I would recommend using the OS Zoomstack which is a nicely styled vector data set that you can reproject and then tile so the pixels don't get "squashed".


-1

I've solved the problem using this: ./configure --with-proj=/usr/local


3

How about something like this: radius=5000 # Maximum distance X=200000 # X coordinate of lookout Y=500000 # Y coordinate of lookout r.import dem.tif output=elev g.region -a rast=elev r.viewshed input=elev output=vshed max_distance=$radius coordinates=${X},${Y} r.to.vect input=vshed output=vshed type=area This should leave you with a raster named "...


0

Got this error too when all my GIS related packages (geopandas, gdal, fiona) that all seem to rely on that libfontconfig randomnly broke. In terminal, I had to first install the fontconfig package using brew brew install fontconfig Navigate to your python install path, which is in your error: cd /Users/tomkom/anaconda3/lib/ Rename the libfontconfig to a ...


2

Here are the replies by mankoff and mmetz from trac.osgeo: I think the issue is related to the region. There is a work-around for the bug: g.region raster=<raster> -a The r.slope.aspect documentation in the NOTES section implies that the region is adjusted to the raster. I do not find that to be the case. Specifically in the NC test data, I've ...


1

Solved as follow: First step: Import HGT data to a temp table: gdal_contour -f PostgreSQL -a elevation -nln contour_lines -i 10 srtm/S01W069.hgt "PG:host=volcano-db user=postgres password=guesswhat dbname=contour" -lco OVERWRITE=YES Step two: Try to COPY the data to the final table. Will fail if we are in the first loop iteration since the table doesn't ...


1

I think the best option for this is going to be using the PGDump driver instead of the PostgreSQL driver, then for the append operations you can use the CREATE_TABLE=OFF layer creation option. The initial command would be something like gdal_contour -f PGDump -a elevation -i 10 srtm/S01W069.hgt | psql --host=volcano-db --user=postgres contour" Then the ...


2

In other words, you want to create a World file from the coordinates of the 4 corners and the width and height of the image 1) You get the width and height of the image with osgeo.gdal, rasterio or any other libraries to open image files as Pillow and others. dataset = rasterio.open('satel.tif') rasterx = dataset.width rastery = dataset.height 2) you ...


0

If we look at the expression you have, A*(A<0.10) (A<0.10) evaluates to a boolean, either True or False If a boolean is used in place of a numeric value, True is equivalent to 1 while False is equivalent to 0 So when A is less than 0.1, we'd expect those values to be multiplied by 1 and remain, while values greater than that would be multiplied by 0....


0

Maybe you use the Chinese characters, I had this problem before. And after deleting the file with Chinese characters, it run successful.


0

I downloaded sample data and took one .dat/.Hdr pair from the archive. Then I had a try with GDAL 3.1.0dev on Windows. This version from the gisinternals comes with a specific driver gdalinfo --formats |find "SNODAS" SNODAS -raster- (rov): Snow Data Assimilation System There is some documentation about the driver at https://gdal.org/drivers/raster/snodas....


0

I solved, @Ian Turton answer helped me to solve the problem. The problem was in two versions of gdal. In the terminal, the environment is set to a newer version of gdal 2.1.3 (/root/miniconda2/bin). CRON runs an older version of gdal 1.10.1 (/usr/bin). I solved this problem by setting up a path to the new version of gdal in the bash script. #!/bin/bash /...


0

An alternative to the approach suggested in the other answers is to use the rasterio package. I had issues generating these using gdal and found this site to be useful. Assuming you have another tif file(other_file.tif) and a numpy array (numpy_array) that has the same resolution and extent as this file, this is the approach that worked for me: import ...


1

The following approach worked pretty well. First I build virtual raster. gdalbuildvrt raster.vrt -srcnodata 0 -input_file_list paths.txt paths.txt is file with following content: a.tif b.tif Then I add a pixel function to it, as showed here https://lists.osgeo.org/pipermail/gdal-dev/2016-September/045134.html. Pixel function is written using numpy, ...


0

The solution to this was quite simple - I was able to open the .gdb file in QGIS and use the 'Select by Expression' function to select the fields I wanted. The formula I used was based on this post: SCINAME IN ( 'kent', 'essex', 'sussex', 'london' ) I was then able to save the selection as a shapefile.


0

I was able to solve this with the help of https://stackoverflow.com/questions/32112597/how-to-crop-polygon-with-latitude-and-longitude-coordinates-from-geotiff-image-i. Ultimately I used the following gdal command line utilities: Create crop shape GeoJson file. I used Mapbox Studio to create my GeoJson. I could have done it by hand, but Mapbox allowed me to ...


0

I am working on a similar problem (wanting to work with rasterized netcdf data) and I think you may have made the same wrong assumption as me, but without information on what instrument your data comes from I can't confirm that. In my case, I assumed the data was already in some kind of grid projection because it came in rows and columns. Instead the 'pixel ...


0

As @Daniel mentioned, you have to specify a value for the masked positions in your masked array. A straightforward way of doing this is calling the .filled() method on said masked array. In your case you could write: # assume arr is the masked array arr = arr.filled(-999) # you can pass an arbitrary scalar to this method.


2

To read an image using gdal as a numpy array you can use the following two lines: ds = gdal.Open('path_to_img', 0) arr = ds.ReadAsArray() However, you can run into a MemoryError if your image is too large. Depending on your needs, you can read one chunk of the image at the time. The ReadAsArray() method can take four useful arguments to accomplish this: ...


1

In your method, change it like the following lines and that will solve the problem: def array2raster(newRasterfn, rasterOrigin, pixelWidth, pixelHeight, array, path): """ This function works fine, until the raster mask. It does the work whatsoever. Not to be used unless necessary. """ # Create a local copy of the array temp = array....


1

Yes there is. you can call the same using GDAL Java bindings. Official documentation can be found at this link For downloading it using maven you can find maven repo link here You can find details of translate function here You can use translate function in two ways public static Dataset Translate​(java.lang.String dest, ...


1

I tried this and have a couple of additional observations. If I drag a geopackage into a qgis project, the raster layers are not named, as was the experience of the OP. However, my geopackage also contains vector layers. These are not imported at all when the geopackage is dragged into the project. However, if I navigate to the geopackage in the qgis ...


3

The code is in GitHub. The Fix Geometries tool is a wrapper to the function LWGEOM_GEOS_makeValid, which then calls different functions for different geometry types. You can read the various LWGEOM_GEOS_makeValid* functions such as the polygon one in the file qgsgeometrymakevalid.cpp


1

There are probably self-intersections in resulting polygon. In QGIS 3 there is Fix Geometries tool. It usually helps running layer through it when there are invalid geometries reported.


1

You need to tell gdal to ignore SSL certificate errors by running the following in the Python Console from within QGIS (from Plugins>>Python Console): import gdal gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES') Unfortunately it looks like you're trying to use an ImageServer, which QGIS isn't really natively compatible with right now and I haven't had ...


1

ArcGIS Pro 2.3 release includes a GDAL/OGR installation. Check all the community ideas that were included in that release here: https://community.esri.com/community/arcgis-ideas/blog/2019/01/17/ideas-in-arcgis-pro-23


3

The easiest way I have found to have GDAL installed on the same machine as ArcGIS is to have a completely separate version of python installed. On this completely separate version of python you should install numpy (pip install numpy) and then download the correct gdal .whl file from https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal This can also be ...


0

You can try var vrt = Gdal.wrapper_GDALBuildVRT_names(vrtFile, aeroFiles.ToArray(), vrtOptions, null, null); vrt.FlushCache(); - it might can help you!


Top 50 recent answers are included