New answers tagged

0

A bit older GDAL versions required the use of -of parameter (notice: with the format name, for example -of PNG) or otherwise it used GeoTIFF as the default format. You can check your version with gdalinfo --version. However, you have given .sid as the extension of the output file and because MrSID driver is read only then any GDAL version will use GeoTIFF ...


1

WarpOptions has xRes and yRes arguments (output x & y resolution). You can either create a WarpOptions object and pass it to gdal.Warp with the options=your_warp_options argument or you can just pass xRes & yRes directly to gdal.Warp as keyword arguments: options = gdal.WarpOptions(xRes=6.5, yRes=6.5, your other args, eg. resampleAlg=gdal....


1

In cases like this I always find that the documentation of the GDAL utility helps me work out what the Python is expecting. So looking at the gdalwarp page leads me to suspect that dstSRS is the same as the -t_srs argument so it's looking for the target SRS (or CRS) - in your case probably the same as the input one.


2

How about something like this, try: gdf = geopandas.read_file(shape) except DataIOError as err: # log error if you want # log(str(err)) raise TypeError clip(raster, shape) Good all try/except Edit: Just to add some links fiona open doc, geopandas read_file is based on fiona open geopandas io doc.


1

I did some steps and with the help of @nmtoken I was able to fix it. I believe that my raster table had some problems maybe because I created it with complex geometries or maybe another reason. This time I changed the raster table I just simply used this codes in PostGIS to create a raster table from a single polygon: SELECT ST_AsRaster(geom, 100, 100 , '...


2

Step 1: Go to https://search.earthdata.nasa.gov/search?q=SPL2SMAP_S Step 2: Select the tile that pops up Step 3: Filter granules by date and region Step 4: Select + Add to Project on right Step 5: Select My Project button that appears Step 6: Select More Options Step 7: Select Customize under Select Data Access Method Step 8: Under Reformat Output select ...


3

In Python, this would be: import gdal import osr driver = gdal.GetDriverByName('GTiff') spatref = osr.SpatialReference() spatref.ImportFromEPSG(27700) wkt = spatref.ExportToWkt() outfn = '/path/to/out.tif' nbands = 1 nodata = 255 xres = 5 yres = -5 xmin = 0 xmax = 680000 ymin = 0 ymax = 1240000 dtype = gdal.GDT_Int16 xsize = abs(int((xmax - xmin) / ...


1

I couldn't open it using /vsigzip//vsicurl/ or /vsitar//vsicurl/. It's a tar.gz of multiple files. Your best option is to download it (manually or you can script with curl or wget, or in python via the requests library). Once you have downloaded it, you can use /vsitar/ to read it: gdalinfo /vsitar/LC08_L1TP_146039_20180119_20180206_01_T1.tar.gz ERROR ...


0

Raster data stored in File Geodatabases is currently not supported by GDAL/QGIS. There is work being done by Nyall Dawson and North Road along with members of the GDAL development community to hopefully bring this feature to light in the upcoming time. Watch this page for news regarding the development of the raster driver: https://north-road.com/blog/


0

With raster RGB color table support : import numpy as np from osgeo import gdal path_inDs = "/data/OCS_2016.extract.tif" path_outDs = "/data/OCS_2016.postpython.tif" driver = gdal.GetDriverByName('GTiff') file = gdal.Open(path_inDs) if file is None: print ('Could not open image file') sys.exit(1) band = file.GetRasterBand(1) lista = band.ReadAsArray(...


0

You need to use gdal.GPI_RGB or gdal.GDT_UInt16 and not gdal.GDT_Int16. That's the idea behind this error !


1

Another option would be to use a Personal Package Archive (PPA). The NextGIS ubuntu PPA currently has GDAL 2.4.0 for xenial, trusty and bionic ubuntu. The ubuntugis stable PPA also has 2.4.0 for bionic currently (January 2020). To upgrade your GDAL to 2.4.0 using the NextGIS PPA, it should look something like this: apt-get install -y software-properties-...


0

I'm not sure, but I think that is not possible. The Attributes Form configurations are stored in the Style of the layer. You can store the style in the same GeoPackage (e.g., saving as Default in the Datasource Database), and the style will be written in the layer_styles table. But styles are XML documents, and a Value Relation widget, in QML format, ...


0

Issue is sorted after creating a copy of the original shape file and running the conversion against the copied shape file. C# consol app >> DataSource shpDatasourceCreateCopy = shpDriver.CopyDataSource(shpDatasource, dataDir + "\CopyTest", new string[0]);


0

This error is extremely confusing, but it is thrown when the in- or output projection is not set. Check your crs and crsGeo to check if they are not an empty. If so, try updating pyproj, setting your PROJ_LIB correctly etc.


1

So here is a full outline of my solution. It seems that the xyz. tiles have an odd formatting since the coordinates are not in the expected order and have a positive N-S pixel resolution. To fix the coordinates I first sorted the coordinates: for file in *.xyz; # sort coordinates do sort -k2 -n -k1 "$file" -o "$file"; done The positive N-S pixel ...


0

I had the same issues with this data. My approach was to create raster-tiles from the point data and then to mosaic it together. However, I did it in R, and I am pretty sure there is a nicer solution, however it worked for me: library(devtools) if (!require(parallelsugar)) install_github('nathanvan/parallelsugar') library(parallelsugar) ...


1

As you can see in following image, I have loaded in QGIS 3 a raster (world_0_360.tif) with longitudes from 0 360 and with latitudes from -90 90. Afterward, I tried out an equivalent command used in accepted answer in your link: gdalwarp -t_srs WGS84 world_0_360.tif world_180.tif -wo SOURCE_EXTRA=1000 --config CENTER_LONG 0 and it ran without any error: ...


1

With the builtin GDAL utilities you can first make the VRT with gdalbuildvrt: gdalbuildvrt mosaic_file.vrt *.tif And then convert to Tiff and reproject with gdalwarp: gdalwarp -t_srs EPSG:32734 -of GTiff mosaic_file.vrt output.tif


1

There is no such create option in GDAL tools but a proper solution would require writing a few lines of Python code. Best workaround that comes into my mind is to select data from any existing data source with point geometries by using SQL that finds nothing. ogr2ogr -of "ESRI Shapefile" -sql "select geometry from test where 1=2" emptytest.shp test.shp ...


1

I would recommend a combination of geopandas and geocube. Here is some untested set of code that should get you pretty close to what you want to do. Step 1: Combine the shapefiles import pandas import geopandas gpd1 = geopandas.read_file("Shapefiles/shp1.shp") gpd2 = geopandas.read_file("Shapefiles/shp2.shp") gpd3 = geopandas.read_file("Shapefiles/shp3....


0

I found my solution. The issue was that I was setting the band.fill(NO_DATA) only when creating the final output. In order to fix this, I simply had to set fill value at the time of rasterizing the layer


0

If you want to keep the 10 largest polygons, I'd recommend executing a SQL sentence on your original data source and then copying the result to a new data source. import ogr ds = ogr.Open('some_file.shp') lyr = ds.GetLayer() name = lyr.GetName() n = 10 # features to keep sql = f'SELECT * FROM {name} ORDER BY OGR_GEOM_AREA DESC LIMIT {n}' # this query ...


3

Assuming your three rasters have the same dimensions, you can use numpy's boolean indexing to accomplish this. First, you need to create three masks, each one corresponding to one of your conditions: con1 = (dcl_array == 1) # raster a is 1 con2 = (dcl_array == 0) # raster a is 0 con3 = (tcd_array == 0) # raster c is 0 Then, you just have to index the ...


0

In the menus above go to raster-->miscellaneous-->build virtual raster This will create the desired merge you want. Then save it to a "real" file.


1

It's my first time writing here, but hope I might give you some insight. I used the code from your post, with some minor changes it worked… kind of. I checked the correlation matrices that I get by comparing some of my classification data with data rasterized in the QGIS tool and raster generated by this code. Raster generated by your code had kappa ...


0

Two options for me : You can do it directly on your raster using connectedComponentsWithStats of opencv and CC_STAT_AREA option. See this link (here). You can do it on your shapefile. First, you have to calculate a new field "area" for your shapefile, then, you can do an attribute query and delete features under a specific value (area = 3 pixels in your ...


1

You need to have GDAL installed. Use before_script call to install GDAL: default: before_script: - apt update && apt install libgdal-dev


0

If you open the tiffs as numpys you can apply any condition you want to them with numpy where. ds = gdal.Open(A) band = ds.GetRasterBand(1) arrA = band.ReadAsArray() ds = gdal.Open(B) band = ds.GetRasterBand(1) arrB = band.ReadAsArray() arrC = np.where((arrA != <yourwhitevalue>), arrA, arrB) [cols, rows] = array.shape driver = gdal.GetDriverByName("...


0

I had the same question when using the gdal_pansharpen.py script. I believed GDAL developped their own algorithm... But this post (here) answered your question... GDAL pan sharpening algorithm = weighted Brovey algorithm Also specified on Github, line 49 (here).


2

In addition to expressions you can use numpy functions in gdal_calc. Sounds like you might be looking for np.where(). Something like: --calc="where((A-B)>0), A, B)"


1

You can create an osr.SpatialReference() object and then use it to set the raster projection using gdal. If you have gdal installed you'll have osr already. Here is an example: from osgeo import gdal from osgeo import osr # create spatial reference object using your EPSG sr = osr.SpatialReference() sr.ImportFromEPSG(32631) # open your raster and set ...


4

As documented in https://pypi.org/project/GDAL/ plain "gdal" is deprecated and you should not use that for writing any new code Additionally, there are five compatibility modules that are included but provide notices to state that they are deprecated and will be going away. If you are using GDAL 1.7 bindings, you should update your imports to ...


0

Maybe your issue is already resolved but while installing QGIS 3.4 of the KyangChaos website I also got stuck on the "known issue" part. On the website it states you have to prepend variable PATH with value: /Library/Frameworks/Python.framework/Versions/3.6/bin:/Library/Frameworks/GDAL.framework/Versions/2.2/Programs: The above doesn't refer to the ...


1

Could you try something like this for python file? C:\PROGRA~1\QGIS3~1.4\apps\Python37\python.exe c:/Users/user/Desktop/xx/run.py


3

You need to set the environment variables to run the QGIS python outside. In this case my QGIS installation is in D:\QGIS, so you should only change OSGEO4W_ROOT. @ECHO OFF set OSGEO4W_ROOT=D:\QGIS call "%OSGEO4W_ROOT%\bin\o4w_env.bat" call "%OSGEO4W_ROOT%\bin\qt5_env.bat" call "%OSGEO4W_ROOT%\bin\py3_env.bat" path %OSGEO4W_ROOT%\apps\qgis\bin;%PATH% ...


10

You are overwriting dst.gpkg every time you run ogr2ogr. From the geopackage documentation: For adding new layers into existing geopackage run ogr2ogr with -update. So do: ogr2ogr -f GPKG dst.gpkg src.shp -nln layerOne ogr2ogr -f GPKG -update dst.gpkg src.shp -nln layerTwo Or if all source shapefiles are in a single directory: Translation of a ...


1

The issue was resolved by the GDAL ticket https://github.com/OSGeo/gdal/issues/2165. The problem was indeed in the type thing. The BigTIFF file has not been created as the standard requires because StripByteCounts is using signed datatype StripByteCounts (279) SLONG8 (17) 98820<21960 21960... whereas acoording to https://www.awaresystems.be/imaging/tiff/...


0

Misunderstanding on my part about DN's interest ... DN=0 (Land), DN=1(Sea) with this command it's ok ogr2ogr output.shp input.shp -where "DN=0" INFO: Open of `../test/test.shp' using driver `ESRI Shapefile' successful. Layer name: test Metadata: DBF_DATE_LAST_UPDATE=2020-01-11 Geometry: Polygon Feature Count: 1 Extent: (-180.000000, 0.000000) - (-...


2

after some research i did the following: gdalwarp -s_srs "+proj=eqc +R=1737400" -t_srs "+proj=ortho +lat=90 +lon_0=0 +R=1737400" input.tif output.tif So basically what caused the error was the -s_srs definition. Since it was the moons dem I had to specify the Radius by +R=1737400. Hope this helps!


5

GDAL supports two SQL dialects which can be used for all data sources. In addition to those two for some databases (Oracle, PostgeSQL etc.) the native SQL can be utilized. The OGR SQL dialect https://gdal.org/user/ogr_sql_dialect.html is made by the GDAL project and it is not at all a full featured SQL language. The other SQL dialect https://gdal.org/user/...


1

gdal_retile -s_srs option sets the source spatial reference system, from your question it is already set as something different. -s_srs is for situations when you know the projection of the data but the data (or associated metadata files) doesn't know or is incorrect. So never use this option if there is a projection set. If you would like to reproject the ...


1

While I'm not familiar with the C#/.NET bindings for GDAL, the library usually divides raster and vector functionality in different namespaces: GDAL (Geospatial Data Abstraction Library) for raster data OGR (OpenGIS Simple Feature Reference) for vector data I assume you need to use the appropriate Open function within the OGR namespace, as you confirmed in ...


0

You may be interested in rioxarray to do this type of operation. It uses rasterio which is an alternative python wrapper for GDAL. import rioxarray import json # load in the geojson file with open("img.geojson") as igj: data = json.load(igj) # if GDAL 3+ crs = data["crs"]["properties"]["name"] # crs = "EPSG:4326" # if GDAL 2 geoms = [feat["geometry"] ...


1

Using the data from ?st_zm, here's a thing with Z and M coordinates: > a Simple feature collection with 2 features and 1 field geometry type: LINESTRING dimension: **XYZM** bbox: xmin: 1 ymin: 3 xmax: 8 ymax: 16 epsg (SRID): NA proj4string: NA a geom 1 1 LINESTRING ZM (1 9 17 25, 2... 2 2 LINESTRING ZM (1 ...


4

This way: >>> from osgeo import osr >>> osr.GetPROJVersionMajor() 6 >>> osr.GetPROJVersionMinor() 2


1

If you are interested in masking the data, I would recommend rioxarray. An example of doing so can be found here. Here is a targeted example for what you probably want to do: from shapely.geometry import mapping from shapely.wkt import loads import rioxarray geom = mapping(loads('POLYGON ((5.5937209 52.24012314 0,5.5936411 52.24048512 0,5.59413417 52....


2

I would do something like this: co = [ 'COMPRESS=JPEG', 'JPEG_QUALITY=90', 'PHOTOMETRIC=YCBCR', 'TILED=YES', 'BLOCKXSIZE=128', 'BLOCKYSIZE=128' ] compressed_ds = gdal.Translate('compressed.tif', 'uncompressed.tif', creationOptions=co) overviews = [2,4,8,16,32,64,128,256,512,1024] compressed_ds.BuildOverviews('GAUSS', overviews) ...


1

Just modify the extent to so that xmin = xmin - 1/2 pixel, xmax = xmax + 1/2 pixel, ymin = min - 1/2 pixel, ymax = ymax + 1/2 pixel. i.e assuming your desired pixel size is 1m and your point layer extent is xmin, xmax, ymin, ymax = 32346264.0,32346299.0,5631990.0,5631999.0 Then set your extent to 32346263.5,32346299.5,5631989.5,5631999.5


4

Another way to do this is to use the rasterio.transform.rowcol() method described in the rasterio transform docs. Example: import numpy as np import rasterio xs = np.array([130.5, 146.0]) ys = np.array([-25.5, -42.0]) with rasterio.open("my.tif") as src: rows, cols = rasterio.transform.rowcol(src.transform, xs, ys)


Top 50 recent answers are included