New answers tagged python
0
To determine hotspots and coldspots, you first need to see globally if there exists some form of autocorrelation/clustering (Global Moran's I).
Then, you could do local statistics to determine local clusters like you did above.
To determine hotspots: Low p-values (significant) areas with Positive z-scores.
To determine coldspots: Low p-values (significant)...
1
To read your raster into a numpy Array you can try:
import gdal
ds = gdal.Open(r'C:\path\to\raster\dsm.tif')
band1 = ds.GetRasterBand(1).ReadAsArray()
2
It can be achieved with a simple walk or listdir:
import os, sys # the standard imports
BaseFolder = r'c:\your\folder\with\data' # change this to match your data
for FullPath, dirs, files in os.walk(BaseFolder):
for ThisFile in files: # iterate the files
fN,fE = os.path.splitext(ThisFile) # separate file name and extension
...
0
May be a bit late but... better late than never right?
this is a good answer i found:
Creating polygon grid using Geopandas
2
You can use gdalwarp to achieve this.
from osgeo import gdal
# open reference file and get resolution
referenceFile = "Path to reference file"
reference = gdal.Open(referenceFile, 0) # this opens the file in only reading mode
referenceTrans = reference.GetGeoTransform
x_res = inputTrans[1]
y_res = -inputTrans[5] # make sure this value is positive
# ...
0
I know this is an old post but someone still might find this helpful. Other than cutting the result with the original you can convert the resulting triangles into centroids, select the centroids inside the original polygon with a inner join query (geopandas.sjoin()). Select the triangles where the intersecting original id is equal to the original polygon id (...
1
This is what I was looking for. I found it in the WKT plugin code
current_directory =
os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
self.action = QAction(QIcon(os.path.join(current_directory, "icons",
"LogoSenara.png")),"&RiegoSenara", self.iface.mainWindow())
2
The latest version of geopandas is 0.6.3. In this version, gdf.crs returns a dictionary like {'init': 'epsg:EPSG_CODE'}. Since it is a dictionary, I think that the most appropriate way is to use tools of geopandas defined in geopandas.tools module.
geom_srid_num = gpd.tools.crs.epsg_from_crs(gdf.crs)
print(geom_srid_num) # OUT: 32616 -> int
EDIT: As @...
1
I would recommend a different approach. Consider the following workflow:
Create a square buffer around a point
Create a second square buffer around a point with a scale factor applied
import pandas as pd
import geopandas as gpd
import geoplot as gplt
# The scale factor and buffer distance
BUFFER_DISTANCE = 3
FACTOR = 0.9
# Make up some data (https://...
2
That's not the way to scale from the center.
You need to scale the difference between the corner and the middle point, and move that scaled difference to the middle:
newX0 = (x0 - middleX) * FACTOR + middleX
newY0 = (y0 - middleY) * FACTOR + middleY
newX1 = (x1 - middleX) * FACTOR + middleX
newY1 = (y1 - middleY) * FACTOR + middleY
2
For #1, yes, they are the same.
For #2, the UTM projection is a system of zones setup for the Transverse Mercator projection. See: https://en.m.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
The string you are using to define the UTM projection is a PROJ string an each parameter is defined here: https://proj.org/operations/projections/...
2
If you are willing to pay for preprocessed data, I believe several of the vendors can be queried for preprocessed data in a specified ROI. My Sinergise sentinel-hub free trial has expired, but think something like that was possible, and a QGIS plugin could help. Not sure what sort of Python API, however.
If you want near- but not-quite-free, you can access ...
1
Using Writing numpy array to raster file as suggested by BERA is the easiest way
With one of my examples
Matplotlib 2D:
Matplotlib 3D (with def axisEqual3D(ax) in set matplotlib 3d plot aspect ratio):
ax.plot_surface(xconv, yconv, zconv, rstride=1, cstride=1, cmap='gist_earth',antialiased=True)
ax.view_init(60,-160)
axisEqual3D(ax)
Convert it into an ...
0
Here is a histogram equalization solution that is working better than linear stretch, but some images are still a little washed out, not as saturated as I'd like. I would love to see if anyone else has a solution.
# Get RGB bands from Landsat image
with rasterio.open(landsat_path, 'r') as f:
red, green, blue = f.read(4), f.read(3), f.read(2)
red[red ...
1
A third way is to create two new variable FOOT and THING in your layer properties then you use @FOOT and @THING in your expression and the calculator will use the value you set in the layer variable.
You can edit the variable to another value and the expression will use the new value. These variable could be used anywhere in QGIS and you only have to edit ...
1
(reviving this coming from a more recent linked Q)
The custom function approach is an effective one, but can be hard to read or maintain if you're not a python wizard. Also, custom functions need to reside in a file in your user profile and so are not automatically portable with the project.
Therefore you could also proceed as follows:
Load the CSV ...
1
When you simply read a geometry value from Microsoft SQL Server, you get a blob in the CLR Type Serialization Format, which is different from WKB.
If you have an SQL Server database, you can cast that value back into the geometry type (although it would be a better idea to export the value as WKB or WKT to begin with).
If you do not have SQL Server, you ...
1
Your indentation is off, but I'll assume that's just a copy/paste issue.
Your issue is that 'F:\\prep\\converted\\*.tif' will try to output a file called literally *.tif which won't work. You need to provide an actual file name. I suggest using os.path.join and os.path.basename
Try this:
from osgeo import gdal
import glob
import os
input_path = 'F:\\...
0
From the author's comment:
I have figured this one out, if you remove the line
pipe.set(renderer.clone()) , it will make an exact duplicate of your
DTM and not as a 4-banded rendered image. Hope this is of some use to
anyone trying to batch save raster tiles like i was.
I can also confirm that this is the solution for the issue where single band ...
1
As noted, reducing the amount of data being called is probably a good idea. To reduce the area you can use a polygon or point geometry. A polygon would use the following:
.filterBounds(ee.Geometry.Polygon([insert coordinate list here ]))
Documentation about how to fill that in can be found at https://developers.google.com/maps/documentation/javascript/...
0
I have solved a similar issue using "gaussian_filter".
from scipy.ndimage.filters import gaussian_filter
data3 = gaussian_filter(data3, sigma=.6)
You can try with different values of sigma.
0
Simple to fix, but took me too long to find this out. Do not set the renderer, i.e., remove this line:
pipe.set(rlayer.renderer().clone())
0
If the only item you have is the .shp file, and you can do without the other parts, consider this solution, opening and saving from OpenJump. https://gis.stackexchange.com/a/306228
You may still need a .proj file, depending on your purpose, but may be able to get away with a "standard" one copied from another project and renamed to match your shapefile, if ...
0
In the end I directly modified the VRT's XML. Rather than using XSLT, I used Python's lxml library with a sprinkling of XPath. I chose lxml rather than ElementTree because of the need to insert CDATA sections.
The following code adds a PixelFunction to the first band:
from copy import deepcopy
from lxml import etree
vrt_tree = etree.parse(tempfile)
...
2
The accepted answer does not work in QGIS3.
Now one has to do:
fieldname='id'
layer=iface.activeLayer()
idx=layer.fields().indexFromName(fieldname)
print(layer.maximumValue(idx))
(I am setting fieldname as the first line to make it easier to cut and paste for someone wanting to test it with their layer)
1
I think your request is simply too large for an interactive session in Python or the JS Code Editor. You may need to do this as a batch task - export the result as an asset and when it completes, import it and then explore it.
Additionally, you may want to break up the analysis by region and date range to make the computation and interpretation of the ...
0
I've worked out a way to do this - although it does not involve deleting all duplicate points (including the 'all but one' normally left by duplicate points.
Instead I have:
Created a buffer around the polygon layer and dissolved the result
Changed the buffered polygon layer to lines
Exploded lines
Ran this code:
layer = iface.activeLayer()
temp = ...
0
You have several way, I use QGis 3.5
In the processing toolbox, go to vector general > delete duplicate geometries. Then select your layer and as output you will have a new "cleaned" layer.
You can also use the Topology checker plugin and then you check for duplicates geometries.
In the first method, QGis will do everything automatically and in the second ...
0
The gdal documentation describes the GDT_Byte as an 8 bit unsigned integer (see here).
So the correct gdal constant is the GDT_Byte. In your code it would be:
driver.Create("out.tiff", cols, rows, 1, gdal.GDT_Byte)
0
Perhaps shapely's unary_union will help, where df is your GeoDataFrame:
import pandas as pd
from shapely.ops import unary_union
hole_indx = pd.isna(df.raster_val)
holes = unary_union(df.loc[hole_indx, 'geometry'])
polys = unary_union(df.loc[~hole_indx, 'geometry'])
polys_and_holes = polys.difference(holes)
0
Not a complete solution, but maybe construct a KDTree from the points in the LineString and then use query to find nearest neighbors to points within the interpolated line.
2
Take a look at np.where and maybe do something like:
y_true= np.ravel(img_r)
y_pred= np.ravel(img_pre_r)
y_true = y_true.astype('int')
y_pred = y_pred.astype('int')
indx = np.where((y_true != -9999) & (y_true != 205) & (y_true != 210) & (y_true != 215))
y_true = y_true[indx]
y_pred = y_pred[indx]
Or apply the np.where call to both true and ...
2
You're right. QComboBox objects allow only single selection because they don't have ExtendedSelection option (as QListWidget objects) for activating multiple selection. However, you can also use a QTableWidget object whose ExtendedSelection option is already activated by default. In following code you have an example.
from PyQt5.QtCore import Qt
class ...
1
Check what gdalinfo says about the image. Notice that you did not get an error but just a warning and it says Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel.
The warning means that the image has, for example, 4 samples for each pixel but the Protometric TIFF tag is set to RGB and there is no ExtraSamples tag ...
0
The accepted answer is not available for QGIS3. In QGIS3, I use the code below:
QgsProject.instance().layerTreeRoot().findLayer(lyr.id()).setItemVisibilityChecked(False)
Also, we can toggle on and off all the layers by:
bool = True # or False
root = QgsProject.instance().layerTreeRoot()
allLayers = root.layerOrder()
for layer in allLayers:
root....
0
When using pyproj, note the differences from various releases in how it is used to transform data. Here are a few examples using new/old capabilities based on the question:
Using pyproj >= 2.2.0
import pyproj
print(pyproj.__version__) # 2.4.1
print(pyproj.proj_version_str) # 6.2.1
proj = pyproj.Transformer.from_crs(3857, 4326, always_xy=True)
x1, y1 = (-...
0
This problem is also caused by changing the version of ArcGIS desktop. I was using ArcMap 10.6 and I had no issues with exporting the jpg file. Recently, I upgraded my ArcGIS to 10.7. Now, its not that the python code does not work. But, every time I open the mxd file try to modify it (even the slightest, like, saving it without even changing anything), I ...
2
One option would be to create a python action for each layer filter you want to set. The action can call the setSubsetString()method on the layer setting the attribute expression. See this q/a for a basic example using setSubsetString():
Subsetting a shapefile and saving it using PyQgis
See link below for how to create a python action.
Creating Custom ...
1
I've looked at that post and It seems very customizable but is quite complex too and It need some windows tools (to execute a batch file through time etc). This is why I propose you other approach just using the QGIS console and a python script.
Here you have a complete script with an scheduled proces. This does a buffer each 10 seconds between an starting ...
1
Solved by looking at this answer:
out_vrt = os.path.join(out_dir, 'mosaic.vrt')
ds = gdal.BuildVRT(out_vrt, [raster1.tif, raster2.tif])
ds = None # this did the trick!
ds = gdal.Open(out_vrt, 0) # 0 = read-only, 1 = read-write.
factors = [128, 256, 512]
gdal.SetConfigOption('COMPRESS_OVERVIEW', 'DEFLATE')
ds.BuildOverviews("AVERAGE", factors)
Explanation ...
2
I agree with @mikewatt, that a good way to handle the error would be to capture it by using try and catch.
But if you want to check that condition, you could do something like this,
def correctPolygon(s):
geoj = s.__geo_interface__
coords = geoj['coordinates']
# check if it is a polygon
if geoj['type'] != 'Polygon':
return False
...
1
Read https://trac.osgeo.org/gdal/wiki/rfc39_ogr_layer_algebra and https://trac.osgeo.org/gdal/wiki/LayerAlgebra and download the Python script from
https://raw.githubusercontent.com/OSGeo/gdal/master/gdal/swig/python/samples/ogr_layer_algebra.py
Running the script without arguments prints help
Usage: ogr_layer_algebra.py Union|Intersection|SymDifference|...
0
This solution worked for me. Just adjust the buffer amount to work with your units.
Solution copied here:
polygon.buffer(10, join_style=1).buffer(-10.0, join_style=1)
2
Using the https://overpass-turbo.eu/ Wizard and some adjustments, I managed to create a query for retrieving the metadata of a specific way (i.e. id:421427136):
/*
This has been generated by the overpass-turbo wizard.
The original search was:
“id:421427136”
*/
[out:json][timeout:25];
// gather results
(
// query part for: “id:421427136”
way(421427136);
)...
1
I am not certain of what you try to achieve but I have two sugestions:
Follow this idea on how to resize a numpy array using scipy.
Use gdal.Warp to warp a tiff to a different resolution. For example:
import gdal
infn = '/path/to/source.tif'
outfn = '/path/to/target.tif'
xres=6.5
yres=6.5
resample_alg = 'near'
ds = gdal.Warp(outfn, infn, warpoptions=...
1
Okay, I figured it out myself. In case somebody cares:
I moved it all to a class, changing
qgs = QgsApplication([], False)
qgs.initQgis()
to
self.qgs = QgsApplication([], False)
self.qgs.initQgis()
Seems to me like the problem was each function was moved to an own thread or something, where it couldn't find the qgs object anymore.
I do have some new ...
1
I found a Python Library called utm 0.5.0.
The conversion.py module in that library seems to do the trick.
I've modified the code by removing the reference to numpy and the custom OutOfRangeError library.
import math as mathlib
__all__ = ['to_latlon', 'from_latlon']
K0 = 0.9996
E = 0.00669438
E2 = E * E
E3 = E2 * E
E_P2 = E / (1.0 - E)
SQRT_E = ...
2
Fiona does not have spatial operations such as clipping, intersections, unions, etc.
The user manual, "explains how to use the Fiona package for reading and writing geospatial data files. "
It states:
There are no layers, no cursors, no geometric operations, no
transformations between coordinate systems, no remote method calls;
all these concerns ...
10
It's not possible to convert GEDI .h5 file to LAS file as including all data. Because that .h5 file includes a lot of information about a point (actually it is a window in GEDI .h5 format, not point). You cannot add all information to LAS file since LAS file has certain attributes for a point not matching attributes/values in that .h5 file. For example, ...
0
Following the answer given by @sbphd, this is what I coded.
import geopandas as gpd
# with columns "id", "latitude", "longitude" - 10k records
df
gdf = gpd.GeoDataFrame(
df,
geometry=gpd.points_from_xy(
df["longitude"],
df["latitude"],
),
crs={"init":"EPSG:4326"},
)
# 10 records
filtered_df
filtered_gdf = gpd.GeoDataFrame(...
Top 50 recent answers are included
