New answers tagged python
0
eventually this worked for me:
dataset = rasterio.open(input_file_path, 'r')
bands = [1, 2, 3]
data = dataset.read(bands)
transform = rasterio.transform.from_bounds(west, south, east, north, data.shape[1], data.shape[2])
crs = {'init': 'epsg:3006'}
with rasterio.open(output_file_path, 'w', driver='GTiff',
width=data.shape[1], height=data....
0
Thank to Nick's answer i was able to get a working code for drawing multiples lines from a coordinates set, not sure if it is efficient but it works. Greetings, Marcelo.
coords = [(0, 0), (0, 10), (10, 0)]
def shp_out(coords):
# IMPORTAR LIBRERIAS
import ogr
import osgeo.osr as osr
# DRIVERNAME
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver....
1
Check out gdalwarp that is part of the GDAL library. There are Python bindings available through the osgeo package.
This is untested, but it should get you close to what you need:
from osgeo import gdal
gdal.Warp(
dstDS,
srcDS,
format='GTiff',
dstSRS='EPSG:3006',
outputBounds=[minX, minY, maxX, maxY]
)
You'll need to define the source and ...
2
You need to iterate through the edges of the LinearRing of the polygon
from shapely.wkt import loads
lin = loads('LineString (289.63171806167395061 -200.22555066079294761, 380.69030837004402201 -65.28898678414094547)')
pol = loads('Polygon ((112.23259911894263041 -229.94933920704846742, 178.75726872246687549 -113.4132158590308137, 309.44757709251092592 -114....
3
You can do this with numpy and rasterio:
import numpy as np
import rasterio as rio
from rasterio import Affine
out_raster = 'path/to/output.tif'
xmin = 147.5
ymax = -34.5
res = 0.00833
cols = 100
rows = 150
profile = {
'driver': 'GTiff',
'dtype': 'uint32',
'width': cols,
'height': rows,
'count': 1,
'crs': 'epsg:4326',
'...
0
For each pair of coordinates:
Create geometry object and add each coordinate (see here).
Create a feature object, set the geometry object and any attribute info, add it to the layer (similar to here)
1
There is another way:
sf = shapefile.Reader(shp=shpFile, encoding='windows-1252')
Probably you'll need to change encoding to an appropriate one.
2
The acceptable types for this type of transaction are listed in the manual as application/xml and application/json not text/xml as you are using.
1
I appreciate your code because it contributes to geodesy, by allowing to make some interesting analyzes in a practical way, regarding spherical approximations for the calculation of length of geodesic lines.
The definition of the spherical radius that approximates the ellipsoid at one point, is a definition that is, at least, ambiguous. What I propose is ...
1
The WMTS server is providing you with a picture of your data, most likely as a png or jpeg file. So, it is unlikely that there are 3 bands or that the values are stored as unsigned ints.
There is no point in looking at the pixel values as they will have been modified by the rendering and compression process. If you want to examine the values of remote ...
2
You could use Result.GetMessage/GetMessages, instead of arcpy.GetMessage/GetMessages
For example:
import arcpy
def log(result):
# Last message
print(result.getMessage(result.messageCount-1))
# Or all messages...
# print(result.getMessages())
log(arcpy.SomeTool_toolbox(some_input))
log(arcpy.AnotherTool_toolbox(some_input))...
1
Based on the article Fran Raga linked to (nyalldawson.net/2016/10/speeding-up-your-pyqgis-scripts ), I did the following:
First:
request=QgsFeatureRequest().addOrderBy('Id')
iter=layer.getFeatures(request)
for festure in iter:
value=feature[fieldname]
This made the features come sorted in the order I wanted them. According to the article, this should ...
2
This is rough, but something like this will work for connecting the polygons in a path. Note that I've change your rectangles to squares for the sake of simplicity, but it should still work. If you need it to work with irregular shaped polygons, you'll probably need to find the centroid.
let currentPoint = [0, 0];
const nodes = [
[-3, 0],
[-3, 3],
[-...
1
Simple indeed....
Proper code below:
from osgeo import gdal
gpkg = gdal.OpenEx('path_to_gpk',gdal.OF_OF_UPDATE,allowed_drivers=['GPKG'])
error = gpkg.DeleteLayer(0)
# this code is also possible as an alternative
# doesn't seem to delete everything though, so possibly DeleteLayer() is better
# gpkg.ExecuteSQL('DROP TABLE [layer name]')
gpkg.ExecuteSQL('...
0
Cracked it!
Basically, the issue was assigning the style to a workspace. After changing my SLD upload code, and the layer publishing code to remove reference to a workspace my code now works as intended.
1
I modified the code from gdal.VectorTranslate returns an empty file
from osgeo import gdal
gdal.UseExceptions()
openOptions = ['ADJUST_GEOM_TYPE=NO']
srcDS = gdal.OpenEx('input2dshp.shp',open_options=openOptions)
ds = gdal.VectorTranslate('output3dshp.shp', srcDS=srcDS, format = 'ESRI Shapefile', layerCreationOptions = ['SHPT=ARCM'])
#Dereference and close ...
0
A solution is to get the DEM elevation values of the same points xy coordinates (red lines projections in the figure below) using osgeo.GDAL or rasterio to read the DEM, GeoPandas for the points shapefile and affine to get the elevation value (see Python affine transforms).
With GDAL
from osgeo import gdal
from affine import Affine
import geopandas as gpd
...
2
There was a matrix algebra error in the Haversine formula. I updated the code in the question. I am getting much better agreement between Haversine and geodesic now:
On my actual dataset:
2
In my mind, the composition of (linear) functions is best interpreted in this way:
Little does that matter to the programming language, so I think it's fine like this:
def transformx(lon, lat):
x = lon*(90-abs(lat))/90
y = lat
return x, y
1
I think DBSCAN clustering algorithm is what you are looking for. If you have PostGIS at your fingertips, then you can simply use st_clusterdbscan function as described here. If you really want python, then you can probably use scikit-learn implementation of that algorithm which is described here.
2
Alphashape returns a Shapely polygon, so if you convert the Matplotlib circle into
a Polygon you can take advantage of Shapely's intersection function.
circle = plt.Circle((50,100), 4, color='r', fill=True)
circle_basepath = circle.get_path()
# Matplotlib's get_path on the circle just gives a circle with radius 1 at the origin,
# so circle_basepath needs ...
0
Try referencing the selected layer by name:
def get_statistic(self):
selected_layer = self.dialog.chooseCombo.currentText()
rlayer = QgsProject.instance().mapLayersByName(selected_layer)[0]
provider = rlayer.dataProvider()
0
This is because newer versions of PsychoPy support both Python 2.7 and 3. The future package is used to ensure the code runs on both platforms.
From command prompt run the following:
pip install future
or via conda command prompt:
conda install future
1
builtins module does not exists in Python 2 (it's called __builtin__). So probably you are using the wrong Python interpreter.
Just to be sure try to print sys.version_info to know which version are you using, and sys.executable to get the absolute path to the interpreter.
I have no experience on Windows but probably you can use the Python 3.7 interpreter ...
0
Abstract classes can be used only as a base class of other classes. There are three classes inherited from QgsGpsConnection. See the class diagram in the docs:
So try to instantiate a QgsNmeaConnection or QgsGpsdConnection or QgsQtLocationConnection.
2
This is really more of a IT/ PyQt question, so I will provide an answer in the context of QGIS plugin development, to try to keep on topic.
The reason your mousePressEvent does nothing is that you have defined it inside your main plugin class. The mousePressEvent method is designed to be re-implemented inside a sub class which inherits from one of the ...
0
if you are using PyQgis 2.x
registry = QgsMapLayerRegistry.instance() urlWithParams =
"crs=EPSG:21037&IgnoreGetMapUrl=1&featureCount=10&format=image/png&styles=&layers=parcels&url=http://localhost:8060/geoserver/cadastralWsp/wms%3FCQL_FILTER=BBOX(geom,
579465,125671,664692,167717)" rlayer = QgsRasterLayer (urlWithParams
, "...
0
You might be interested in rioxarray.
import rioxarray
xds = rioxarray.open_rasterio(...)
https://corteva.github.io/rioxarray/html/examples/reproject.html
xds_lonlat = xds.rio.reproject("epsg:4326")
https://corteva.github.io/rioxarray/html/examples/clip_geom.html
clipped = xds.rio.clip(features, features_crs)
0
Esri ArcGIS can read netCDF files directly. No need to convert to ascii.
0
Just in case anyone wanted to see what I did, here are the changes to the code that allowed it to work:
with collection(os.path.join(newShpDir, outputShp + ".shp"), "w", "ESRI Shapefile", schema) as output:
# Merge all of the features into one feature
for i in range(len(features)):
try:
merged = cascaded_union(features[i])
...
1
I use Mac and QGIS 3.4 and got the same error.
After checking the PyQt5 folder in QGIS, I found there is no QtChart file in the folder, which is the reason why we cannot import it.
And the version of PyQt5 in QGIS is 5.11.3
The folder for PyQt5 in QGIS is
/Applications/QGIS3.app/Contents/Resources/python/PyQt5
I have an anaconda version of PyQt5, which ...
3
Your example is a MULTIPOINT.
Members of a multi-point collection are accessed via the geoms property or via the iterator protocol using in or list().
Here's an example:
>> from shapely import wkt
>> from shapely.geometry import MultiPoint
>> p = wkt.loads("MULTIPOINT (92 169, 100 163.552380952381, 105.2167832167832 160, 266....
2
I'm not sure why the order matters to you but the shapely manual shows wrapping the MultiPoint in a list:
pprint.pprint(list(points.geoms))
so I would expect something like:
points = MultiPoint(list(points.geoms)[::-1])
to work. Though the manual makes no claims that MultiPoints honour the order of the input list.
4
Yes, you can write a one bit raster with rasterio*.
You need to:
write to a format that supports a 1bit dataset, such as GeoTIFF;
ensure your numpy array is np.uint8/ubyte so rasterio doesnt raise the TypeError: invalid dtype: 'bool' exception; and
pass the NBITS=1 creation option to tell the underlying GDAL GeoTIFF driver to create a one bit file.
import ...
0
There is some sort of bug in ArcPro. The work around is this:
Make the script output "Desired Name" a precondition to the FeatureClass to FeatureClass tool. Do not connect it to the parameter.
Open the FeatureClass to FeatureClass tool and type into the Output Feature Class Parameter %Desired Name%.
That worked for me.
4
If you call rasterio.dtypes.check_dtype(np.bool_) you'll see that it's not a known dtype, because gdal doesn't support a true 1-bit dtype. GDT_Byte is the smallest. The list that rasterio is checking against is:
dtype_fwd = {
0: None, # GDT_Unknown
1: ubyte, # GDT_Byte
2: uint16, # GDT_UInt16
3: int16, # GDT_Int16
4:...
1
Don't know if this is the simplest way, but it saves gathering elevation data. The USGS-National Map has a REST service that you can use to query elevation for lat/lon coords.
Service url:
https://nationalmap.gov/epqs/
You can use pythons requests library and format your query string according to the service parameters. You need your input coordinates in ...
0
I found an answer here. There are several examples here, but this one works the best for my case:
source = ogr.Open('MyGeoPkg.gpkg')
layer = source.GetLayer()
schema = []
ldefn = layer.GetLayerDefn()
for n in range(ldefn.GetFieldCount()):
fdefn = ldefn.GetFieldDefn(n)
schema.append(fdefn.name)
print(schema)
4
a shapefile can only contain one geometry-type. For each of them, you have to create another shapefile.
More info: https://gis.stackexchange.com/a/74629/7849
2
From the linked article (emphasis mine):
The ArcGIS Server component of ArcGIS Enterprise also ships with Python. In recent releases, both Python 2.x and Python 3.x runtimes are provided.
Below that, it is specified that versions 10.6 and higher come with Python 3.6.x and 2.7.x. I see no indication that 3.6 is not included in any 10.7 version.
ArcGIS ...
1
Yes you can pass an array. The documentation specifies:
rasterio.plot.show(source, etc...)
Parameters
source (array or dataset object opened in 'r' mode or Band
or tuple(dataset, bidx))
Yes, it's the same. Demo using rasterio.plot.adjust_band that show uses to do the adjustment:
import rasterio.plot as rp
import numpy as np
def norm(...
1
You effectively want to iterate over columns of GeoDataFrame, which you can get using geopkg.columns.
So assuming you don't want to iterate over geometry field:
for col in geopkg.columns.drop('geometry'):
data = geopkg[col]
0
I found the answer - I found the .pth file located within the ArcGIS Pro installation itself (default location C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\site-packages\ArcGISPro.pth) and used it as the basis for the .pth file in my virtual environment, which now looks like this:
import os; d = r"C:\Program Files\ArcGIS\Pro\bin"; os....
0
Use the area property of the geometry column:
import geopandas as gpd
df = gpd.read_file('/home/bera/someshapefile.shp')
df['area'] = df.geometry.area #/10000 for hectares
0
For those of you interested in the future, divinding by the quantification value as suggested by @Zac Wang produces the right result
Another option is to change the data type from uint16 to float. Using the following code, the issue is solved:
red = rasterio.open(S_files[2]).read().astype('float')
nir = rasterio.open(S_files[3]).read().astype('float')
ndvi ...
0
In case someone is interested in the future, I managed to make it work with the following code.
First, open an empty list. Then import the shp with geopandas. Next, we access the geometric information and extract it. The output of that is a dict containing the information required by rasterio.mask.mask
geom = []
StudyA = gpd.read_file('.../...
2
Revised Answer: The Path settings are stored in the qgis-bin.env file in the bin folder. If you execute the qgis-bin.exe file directly, bypassing the Path settings in the qgis.bat file, QGIS still has the same Path because it is loaded from the .env file. (I don't know what the interaction is between the .env Path and .bat Path, whether one trumps the ...
0
You seem to be on the right track. On my computer I found the arcgisscripting library here:
C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages
perhaps you could try your .pth file with just two entries:
C:\Program Files\ArcGIS\Pro\Resources\ArcPy
C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages
see ...
0
The error message is pointing out your problem. The geometry filter of Planet Data API does not support FeatureCollection. It only supports 'Point', 'MultiPoint', 'LineString', 'MultiLineString', 'Polygon', 'MultiPolygon'. Here you can check a basic example.
As a workaround, you could loop throught your feature collection and apply the API call for every ...
0
I got a similar weird result when I used complete jp2 images (with DN -digital numbers- values without transforming in reflectances) for calculating NDVI index for a Spanish region.
import rasterio
b8 = '/home/zeito/pyqgis_data/jp2/T30STF_20170422T110651_B08.jp2'
b4 = '/home/zeito/pyqgis_data/jp2/T30STF_20170422T110651_B04.jp2'
NIR = rasterio.open(b8).read(...
Top 50 recent answers are included
