New answers tagged

0

You should be able export RGB JPEG using gdal_translate with [-expand {gray|rgb|rgba}] flag. You should use command like this: gdal_translate -of jpeg -expand rgb -a_nodata 255 -co worldfile=yes col.vrt col.jpg There is included -a_nodata to leave nodata pixels white not black and -co worldfile=yes to include georeferencing file. However I'm not sure about ...


1

Thanks to the extensive help of @mdsumner, this solution works by using the polymer R package to break polygons down into a mesh of triangles, calculating the overlapping triangle segments, and then reassembling the triangles into polygons. This approach is slower but it seems to generally be robust to these non-noded intersection errors which here stem from ...


0

Faster solution with just 1 layer definition: <OGRVRTDataSource> <OGRVRTUnionLayer name="DTM_Points"> <OGRVRTLayer name="T_test"> <SrcDataSource>CSV:T_test.txt</SrcDataSource> <SrcSQL dialect="sqlite">SELECT field_2, field_4, IFNULL(field_7, IFNULL(field_6, field_5)) AS Z ...


2

To combine each of the TIFF files into one multiband TIFF file, you should use gdalbuildvrt. The command would be: gdalbuildvrt -separate -resolution highest -r cubic combined-raster.vrt ...B01 ...B02 ...B03 ...B04 ...B05 ...B06 ...B07 ...B08 ...B8A ...B09 ...B10 ...B11 ...B12 This would provide a single VRT-raster, which can then either be used directly, ...


3

Read the documentation of VRT https://gdal.org/drivers/raster/vrt.html and add color table. Step by step example: Create test image with gdal_create https://gdal.org/programs/gdal_create.html#gdal-create. Utility is included in GDAL version 3.2 and higher. gdal_create -of GTiff -outsize 10 10 -bands 1 -burn 2 -ot Float32 float_gray.tif Create a VRT from ...


0

I made a test with the json that fails for you. I mean this one { "type": "FeatureCollection", "features": [ { "type": "Feature", "geometry": { "type": "Polygon", "coordinates": [ [ [ 182....


0

Okay, I haven't found any suitable solution for splitting custom polygon by the antimeridian. I implemented it by myself using Pyproj and Shapely. Maybe someone will find it useful. from pyproj import CRS, Transformer from shapely.geometry import Point, LineString def intersects_antimeridian_precise(start_point, end_point): ''' Detects if line ...


0

I am using Windows 10 and the .SAFE format seems like a simple folder to me. You can find the *.jp2 raster files in this folder: *\*filename*.SAFE\GRANULE\*filename*\IMG_DATA You can combine different bands (even with different pixel sizes) using the r.composite available in the GRASS toolbox (just tested it with QGIS 3.14).


0

Build virtual raster Right-click and Save/export


1

Following comments from @user2856 and @user30184 and info here, the solution was: gdal_translate -of NetCDF -a_nodata -9999 -a_ullr -17367530.45 7314540.83 17367530.45 -7314540.83 -a_srs "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" infile.h5 tmp.nc gdalwarp -overwrite -s_srs EPSG:6933 -t_srs ...


0

On linux, the epsg database is in your system's proj4 folder (for example /usr/share/proj/epsg) as a simple text file. You can add the appropriate entry in wkt format, providing you have root access.


0

I added this pre-processing data from netCDF to tiff with gdal binary, in case, someone needs: Example OLCI level 2 data: S3A_OL_2_LFR____20180131T225040_20180131T225340_20180202T040253_0180_027_215_2520_LN1_O_NT_002.SEN3 + First, create the lat.vrt file: gdal_translate -of VRT NETCDF:"/home/rasdaman/RASDAMAN_PROJECTS/PARSEC_MUNDI/Sentinel-3/...


1

As @Hornbydd suggested for rivers it is way better to download existing dataset. E.g. from hydrosheds or OSM (waterways and water layers). And clip that dataset to your study area. You can relatively easily extract drainage. E.g. with Grass tools as suggested in other answer. If that is all you need, good. But drainage is far from rivers, depending on detail/...


2

I looked into this and was not able to reproduce the exact results you got from SNAP (specifically I'm not sure why your output has a minimum value of 5). However, I was able to confirm that the results from SNAP and GDAL differ significantly. Two insights: GDAL is implementing the 95% clip correctly. I examined the original raster to identify the minimum/...


3

There are various implementations of watersheds and hydrology delineation using the the A* algorithm. For instance in grass gis r.watershed and python pysheds. In this instance I recommend the grass gis implementation as the method was developed on SRTM data.


2

gdal.Info might be the simplest way: ds = gdal.Open('/vsicurl/https://github.com/OSGeo/gdal/raw/master/autotest/gcore/data/3376.tif') info = gdal.Info(ds, options='-json') print(info.keys()) dict_keys(['description', 'driverShortName', 'driverLongName', 'files', 'size', 'coordinateSystem', 'geoTransform', 'metadata', 'cornerCoordinates', 'wgs84Extent', 'rat'...


2

Indeed, as @snowman2 points out, using pyproj fixes the performance issue. The relevant command would look like this (for more complex geometries use shapely.ops.transform): python -m timeit -s "from pyproj import Transformer" -s "transform = Transformer.from_crs(31287, 4236).transform" "transform(419908, 333400)" It sets up a ...


3

I would recommend using pyproj as it has dealt with this issue already: https://pyproj4.github.io/pyproj/stable/advanced_examples.html#optimize-transformations The creation of the transformer has more overhead in PROJ 6+. That is why pyproj added the Transformer class. See: https://github.com/pyproj4/pyproj/issues/187


0

Here is the python gdal code to extract overview from Sentinel 2 JP2 format with overview level 2 (0-index). import gdal jp2 = "T32UMD_20190727T104029_TCI_10m.jp2" dataset = gdal.Open(jp2) band1 = dataset.GetRasterBand(1); overview1 = band1.GetOverview(2); band2 = dataset.GetRasterBand(2); overview2 = band2.GetOverview(2); band3 = dataset....


0

For this type of work I would get the EGM96 or EGM2008 elevation raster data. Then if you want to go just with GDAL, look at the LocationInfo function to extract for each pair of XY locationsthe corresponding EGM value. https://gdal.org/programs/gdallocationinfo.html I believe you can get both EGM models in raster format from: https://earth-info.nga.mil/...


0

I came across the same issue and created a AUR package for gdal with ecw support: https://aur.archlinux.org/packages/gdal-ecw It's based on the solution above.


1

These kinds of problems are due to the fact that the current GEOS overlay (intersection in this case) algorithm is not totally robust. It can fail on some data, typically that which contains nearly-coincident linework. It is probably that the buffers do contain these situations, if the original points are close. Happily, a much-improved overlay algorithm ...


1

GDAL does not invent the EXIF metadata for you. When you run gdalwarp -t_srs EPSG:4326 IW32001452_mask.tif /tmp/tif2jpeg.tif if IW32001452_mask.tif does not contain EXIF tags then tif2jpeg.tif does not have them either and even -co write_exif_metadata=yesprepares the JPEG file to include EXIF metadata, if there does not exist any tags to write then there ...


1

Got some help and figured it out. Changed some things around and basically need to use gdal.ReadDir('/vsimem/temp/') to get a list of all files created and then add those files one by one to s3.


1

Finally I managed to make a simple polar projection without georeferencing the raster image. For that I used ImageMagick with the following command line: convert raster_input.jpg -virtual-pixel HorizontalTile -background Black -distort Polar 0 -fuzz 50% -trim polar_output.jpg I know there is another solution with Gimp through the Filters > Distort > ...


1

It is not very difficult to do that with python gdal module. You only need two X references for avoiding rounded shapes. First, I downloaded your image and arbitrarily assigned a CRS to it. It looks like (for one band) as follows: Afterward, I put two reference points (in blue) as in following image. I write down its respective X coordinates (first_Xref = ...


0

You may have to replace dstSRS='EPSG:4326' with gdal.WarpOptions(dstSRS='EPSG:4326') to get rid of the error that you are getting. You are supplying a keyword object where a gdal.WarpOptions object is needed.


0

If setting the GDAL_DATA environment variable does not work then you may have to replace dstSRS='EPSG:4326' with gdal.WarpOptions(dstSRS='EPSG:4326') to get rid of the error that you are getting. You are supplying an keyword object where a gdal.WarpOptions object is needed. Check out the docs for more information: https://gdal.org/python/osgeo.gdal-module....


0

Even though it is in beta (on CRAN) I would recommend migrating to the terra package. The terra package is from the same developer as raster and is considered an eventual replacement for raster. The functions/code have mostly been moved over to C++ and exhibit massive performance and speed gains, including in extract. One thing to note is that vector (points/...


0

Finally, I found the answer. Why it cost too long is that the GDAL(3.1.2) used for test does not support spatial index at all. Since the docs said: Since GDAL 3.2, the driver can use the native .spx spatial indices for spatial filtering. Why I did not notice the version was that I thought the 3.2.x would be the current stable version since the docs ...


0

If you don't want to use environment variables (or dont have permissions) adding --config MSSQLSPATIAL_USE_GEOMETRY_COLUMNS NO to the command has solved this issue for me.


0

Have you tried to add the s_srs parameter? Input epsg is 4326 but perhaps the geotiff CRS definition is not the standard expected from GDAL. Moreover are you sure about the t_srs parameter? Are your data on UTM12N zone?


1

@Jarekj1 yes, I also noticed it trying to change a plugin from QGIS2 to 3, in QGIS3.4. I could notice that in QGIS2 it accepted simple geometries, but in 3 all of them are "apparently" generated as multipart. Just in case, I put a conditional in my code if feature.geometry().isMultipart(): to handle the case that a simple geometry arrives, but it ...


2

It is hard to give an objective answer. Generally speaking OpenFileGDB driver's main purpose is to help GIS users make data accessible without need of proprietary software - read only. In other words fix accessibility problem. OpenFileGDB driver was created by reverse-engineering, and so might not / does not have some functionalities as good as native ...


0

As documented in https://gdal.org/programs/gdal_grid.html -txe <xmin> <xmax> Set georeferenced X extents of output file to be created. -tye <ymin> <ymax> Set georeferenced Y extents of output file to be created.


0

As a note for spack/module users: after installing fiona spack install py-fiona You should load both the modules for py-fiona and gdal (if using module load) to avoid getting the ImportError: libgdal.so.26. gdal should also have been automatically installed as a dependency of fiona.


0

That warning doesn't affect your output at all. As the other folks suggested, sf package can be useful for cleaner and simple EPSG coding. Here is your code adapted: library(sf) library(mapview) coords <- list(rbind(c(1,3), c(2,8), c(3,10), c(4,12), c(1,3))) x2 <- st_sf(st_geometry(st_polygon(x = coords)), crs="EPSG:4326") x2$ID <- 1 ...


1

Create a tool that does what you require in ArcMap or ArcGIS Pro. You can create the tool either with Model Builder or using a Python script. The tool should have an OUTPUT PARAMETER for the location of the output file generated. Test it to make sure that it generates the file that you require in ArcMap or ArcGIS Pro. After a successful run of the tool, ...


2

I was also quite confused about the move away from PROJ Strings and what to use instead until I watched this video of a talk held last year at FOSS4G. Here's my main takeaways: Why PROJ and GDAL are moving away from PROJ-strings With PROJ4 every transformation is done using the "hub-approach". So everything will first be converted to WGS84 and only ...


2

I had a look into gdal2tiles.py code, and there is commented out part in update_no_data_values function with TODO: gbataille - check the need for this replacement. Seems to work without replace BandMapping tag for NODATA bands.... Well it seems it is not working that well. For now the fastest workaround I could find: In QGIS top menu use Raster / ...


3

I managed to do it, but I think it's unnecessarily complicated: I used "snap points to points" (saga_cmd shapes_points 18) @Taras suggested, to snap one of my layers to the other one. Then, since I am not aware of SAGA having a spatial join tool, I created small buffers around the points which stayed at their location with saga_cmd shapes_tools ...


Top 50 recent answers are included