New answers tagged python
0
Petar,
here's my proposal :
You open the .csv file ("add delimited text layer" - CRS : WGS 84 - X and Y being the longitude and latitude attributes).
You can save it as a shape file or any format if necessary.
In the Processing Toolbox you choose "split vector layer", as "unique ID Field" you choose "DAY" and the toolbox generates the awaited files ...
...
3
Be careful with your variable names. You're dangerously re-defining the variable c within your list comprehensions.
Do this instead:
counties = [x['long_name'] for x in c if 'neighborhood' in x['types']]
localities = [x['long_name'] for x in c if 'locality' in x['types']]
Better yet, use different names entirely so your code is easier to read and reason ...
0
passing lat / lon directly as a string tuple fixed the problem - didn't even need to convert to shapely point to turn into LineString
geometry = [xy for xy in zip(linkdf.lon, linkdf.lat)]
linkdf['geometry'] = geometry
linkdf_line = linkdf.groupby(['linkId'])['geometry'].apply(lambda x: x.tolist())
linkdf_line = linkdf_line.apply(lambda x: LineString(x))
...
0
SpatialReference.GetAuthorityCode() takes None as a parameter, that finds an authority node on the root element (ie. projected/geographic as appropriate). Same applies to GetAuthorityName():
In [1]: import osgeo.osr as osr
In [2]: srs = osr.SpatialReference()
In [3]: srs.SetFromUserInput("EPSG:27700")
Out[3]: 0
In [4]: srs.GetAuthorityCode(None)
Out[4]: '...
3
Creating a point with 2 coordinates actually creates a 3D point as can be tested with
point=ogr.Geometry(ogr.wkbPoint)
point.AddPoint(float(1.0),float(2.0))
print point.ExportToWkt()
POINT (1 2 0)
You can turn 3D points into 2D points with FlattenTo2D
point.FlattenTo2D()
print point.ExportToWkt()
POINT (1 2)
1
gdal supports the GRIB format. You can use it to open your file and read a specific band as a numpy array. Then, assuming both your grib file and the x y coordinates share the same spatial reference and that your raster is not skewed, you can calculate the row and column (indices) of the array the coordinates correspond to. To do this you need to know the ...
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
If you modify your code to create a list of True/False instead of geometries/None you should be able to use loc and the list as mask:
Access a group of rows and columns by label(s) or a boolean array.
point_nodes = [True if x is not None else False for x in point_nodes]
newships = ships.loc[point_nodes]
Which will give you all columns of ships but only ...
1
The error message says that the key 100 is not found in this line uri = TYPE_MAP[layer.wkbType()]. So layer.wkbType() returns 100 which is the code for NoGeometry. So I guess you have a layer with no geometry set. If you can somehow disable or remove that layer it might work. Alternative try to update QGIS2Web. Looking at the source code you can see that an ...
0
With the solution of How to find which points intersect with a polygon in geopandas? and your data
....
polygon = sh.polygon.Polygon([[p.x, p.y] for p in prt.geometry.tolist()])
polygon.wkt
'POLYGON ((196320.0493430012 -479448.5877129403, 199780.1489409541 -446899.1169124786, 161437.1203666388 -440485.6615666826, 192413.1253873391 -465423.3888242305, 196039....
2
If you need modify the "Feature" Blending Mode on "Darken" use this:
layer.setFeatureBlendMode(QPainter.CompositionMode_Darken)
2
ESA's SNAP software does have a Python API which you can set up by following these instructions.
Additionally, in your SNAP bin directory you'll find the Graph Processing Tool executable. This provides a command line interface to the operators available in the SNAP GUI. Further, you can create your own graphs depending on your workflow and specify your ...
0
Based on your input bounds, it appears that the latlon bounds are likely a mask.
The shape of the latlon bounds are not rectangular in the latlon form:
Or in the projected form:
Also, if you create a transform from the bounds:
The xscale and yscale are ~2.6 and ~1.3 respectively. Which is much smaller than the expected 2k-3k. So, the data is likely a ...
2
Check out PySAL: The Python Spatial Analysis Library
https://pysal.readthedocs.io/en/latest/
0
You actually need to change this on the resource level:
layer=cat.get_layer("layer_name")
resource=layer.resource
resource.projection='EPSG:4326'
cat.save(resource)
cat.save(layer)
cat.reload()
1
I think you are trying to split the polygon with the same polygon. You only use poly.geojson
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 ...
0
I would recommend upgrading to pyproj 2.2.0 as it should take into account the datum shifts.
https://pyproj4.github.io/pyproj/v2.2.0rel/examples.html
>>> import pyproj
>>> pyproj.__version__
'2.2.0'
>>> pyproj.CRS("EPSG:27700")
<Projected CRS: EPSG:27700>
Name: OSGB 1936 / British National Grid
Axis Info [cartesian]:
- E[...
1
The "empty" image contains QA bands that are completely masked out. You can avoid the empty sample by first selecting the bands that have valid data. For example:
var training_empty = image_empty
.select([
'B1',
'B2',
'B3',
'B4',
'B5',
'B6',
'B7',
'B8',
'B8A',
'B9',
'B10',
'B11',
'B12',
// 'QA10',
...
0
The solution we used in the end was to adapt scripts used for creating the yamls for Sentinel/Landsat data which use Gdalinfo to create the yaml files. The updated product yaml is below, although it seems to have little different from the one linked in question.
The difference I suspect is that the product extents are to a far higher degree of precision, ...
0
Duplicate of: https://stackoverflow.com/questions/54972600/get-central-meridian-pyproj
I believe your answer is in pyproj 2.2.0 with the CRS class: https://pyproj4.github.io/pyproj/v2.2.0rel/api/crs.html
>>> import pyproj
>>> pyproj.__version__
'2.2.0'
>>> crs = pyproj.CRS("ESRI:102719")
>>> crs
<Projected CRS: ESRI:...
1
pyproj 2.2.0 can solve this:
>>> import pyproj
>>> pyproj.__version__
'2.2.0'
>>> transformer = pyproj.Transformer.from_crs("esri:102001", "epsg:4326")
>>> transformer.transform((-5171461.906673405, 3372679.8809247203), (209275.49942131247, 5395930.584568907))
((22.775277555142978, 62.48292643585552), (-145.4797312671488, ...
0
I had forgotten to add the Header
#!/usr/bin/python
# -*- coding: utf-8 -*-
before I added the script to the toolbox, I added the header after and it didn't affect the encoding until I deleted and then added the script to the toolbox again,
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
Edit: After reading your subsequent comments above about what you are trying to achieve, I now doubt this actually what you want...
I have written the code below which you should be able to adapt. I have tested this in QGIS 3.4 and have included some commented lines which should allow for adaptation to QGIS 2.18. No guarantees though as I have not used the ...
0
The QgsVectorLayer class does not have the function geometry. https://qgis.org/api/classQgsVectorLayer.html. Maybe you should use:
for feature in layer.getFeatures():
print(feature.geometry().length())
1
Let’s point out some useful information:
1 - A layer has features.
2 - A feature has a geometry
With that in mind, you need to iterate over the features in your layer. For each feature you do:
for feature in layer.getFeatures():
geom = feature.geometry()
len = geom.length()
And so on...
1
You have a problem with your shapefile because the formulation is correct
source = driver.Open("test.shp", 1)
layer = source.GetLayer()
# numbers of records
len(layer)
43
# the fields of the shapefile
layerDef = layer.GetLayerDefn()
print([layerDef.GetFieldDefn(i).GetName() for i in range(layerDef.GetFieldCount())])
['DIP', 'DIP_DIR', 'RHR']
# first ...
2
You don't need to use a cursor for this.
You can create a layer using a 'where' clause so that it only includes the features that match your criteria. Then copy that layer to your new shapefile.
If fact, you can even skip creating a layer and you can use FeatureClassToFeatureClass_conversion() to do it all in one step as that function includes an optional '...
0
Update: I figured out a good way to do this using the Basemap package in Python. I wrote up a tutorial that walks through how to plot several .shp files on different maps all with the same extent. It is posted at https://mapping-with-python.herokuapp.com/index.html
The essential components are to open a map object
map = Basemap2(llcrnrlon= LongLeft,...
1
I hope you have an answer to this by now, but just in case not I would suggest adding a line after the imports that sets the prefix path. It should look something like os.environ['QGIS_PREFIX_PATH'] = r'C:\OSGeo4W64\apps\qgis'. The code here is for Windows, but it should be similar in macOS. This worked for me.
0
You can use a similar method for the marker but by adding symbol.symbolLayer(0) to identify the initial simple marker layer:
layer = iface.activeLayer()
single_symbol_renderer = layer.renderer()
symbol = single_symbol_renderer.symbol()
symbol.symbolLayer(0).setStrokeColor(QColor(255,255,1))
symbol.symbolLayer(0).setStrokeWidth(0.9)
layer.triggerRepaint()
...
0
I am not sure this is the solution you are looking for, but I think it works. If you want to show the change between rasters, this accomplishes that goal. This is the workflow for comparing the change between 2 rasters.
The workflow is:
1. Convert rasters to vector using the GDAL: Polygonize tool
2. Use the QGIS: Difference tool to compare change(found ...
2
This line is essentially overwriting the entry in the ftcs dictionary within the loop:
ftcs[row[0]] = [row[1], segment_dates.get(row[1])]
Hence you ending up with just one date entry.
I think you need to explore the dictionary method has_key() to check for the existing entry then return the value (your list) with ftcs[row[0]] ,lets call it mylist, update ...
2
Apparently, setting table config on QgsAttributeTableView has no effect on column visibility. I have to set the column config on QgsAttributeTableFilterModel. This works for me:
def setHidenColumns(self, *args):
self._config = self._layer.attributeTableConfig()
columns = self._config.columns()
for column in columns:
if column.name in ...
0
I use rasterio and geopandas. My example uses UTM coordinates. Obviously these fields will depend on your particular shapefile. In my experience this produces indentical results to the QGIS Point Sampling Tool. I like this method because the resulting DataFrame of point and corresponding raster values is easy to analyze (e.g. compute the difference between ...
2
If you have a SpatiaLite database you can use the SQL function ST_Is3D http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html. Example that returns "1" as a result
SELECT ST_Is3d(ST_GeomFromText('POINTZ (1 1 1)'));
CoordDimension and ST_NDims can by used as well, the first returns "XYZ" and the latter "3" but I think that ST_Is3D is the best match for ...
2
what you're describing is called ‘string interpolation’.
In days of old, you could simply add strings together with the + operator, but this looks a bit ugly.
As you're using Python 3.6, the easiest way is called “f-strings”, which is new to Python 3.6
I admit I’ve not tried this, but something like this should do the trick. Put an 'f' before the string, ...
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 ...
0
A solution that works for me it's to divide script into two:
- first script iterate on an item list, pass an "id_item" to the second, and erase temporal files
- second script do geoprocess for each item.
Note: in sample code there is a comented sample geoprocess, it's only for take an idea of what could be.
First script (01_call_process.py):
import ...
1
I don't know if this is The Right Way of handling this, but I came up with the following solution:
import io
import fiona
import geopandas as gpd
def read_tracks_filtered(path):
src = fiona.open(path, layer='tracks')
meta = src.meta
meta['driver'] = 'GeoJSON'
with io.BytesIO() as buffer:
with fiona.open(buffer, 'w', **meta) as dst:
...
0
All calls to PostGIS functions must be schema qualified: schema_name.function (source).
To bypass writing the schema every time a PostGIS function is used, map the schema where PostGIS is (probably public) to the search_path (see here). As admin:
ALTER DATABASE <database_name> SET search_path TO schema1,schema2;
Moreover, make sure the private user ...
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
The easiest way is inputting the GeoJSON URL directly into gpd.read(). I'd tried extracting a shapefile from a zip before this using BytesIO & zipfile and had issues with gpd (specifically Fiona) accepting file-like objects.
import geopandas as gpd
import David.SQL_pull_by_placename as sql
import os
os.environ['PROJ_LIB'] = r'C:\Users\littlexsparkee\...
0
Based on @joris' comment, I updated the solution there to current versions:
thisFile = "yourFileRoot"
with rasterio.Env():
with rasterio.open(thisFile+'_AVE_DSM.tif') as src:
with rasterio.open(thisFile+'_AVE_MSK.tif') as msk:
image = src.read(1) ## first and only band
mask = msk.read(1) ## first and only band
...
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.
1
Try this code. It creates a sorted csv (csv_sorted) and appends two new columns: 'time_seconds' and 'trip_id'.
import pandas as pd
import numpy as np
def getSec(time):
return sum(x * int(t) for x, t in zip([1, 60, 3600], reversed(time.split(":"))))
# Read in csv
my_file = r"C:\Users\Jon\Desktop\output2019_05_25-11.csv"
csv = pd.read_csv(my_file, ...
1
You can use collections.defaultdict(list):
from collections import defaultdict
import arcpy
FTC = r'C:\data.gdb\features'
FIBERCABLE = r'C:\data.gdb\cables'
d = defaultdict(list)
with arcpy.da.SearchCursor(FTC, ['conduitipid', 'segmentid']) as cursor:
for con, seg in cursor:
d[con].append(seg)
d2 = {k:v for k,v in arcpy.da.SearchCursor(...
2
Ways are composed of nodes as you surmise. You can get geometry by modifying the out call (see the documentation):
bb gives bounding box
center gives a centroid
geom gives ful geometry (probably the option you want and possibly with the JSON option).
2
The most efficient way is to create two dictionaries:
Map segment ID to placed/planned date
Map conduit ID to your desired list by using the first dictionary
Here's an example:
import arcpy
tbl_cables = 'FIBERCABLE'
tbl_ftc = 'FTC'
field_segment = 'segmentid'
field_conduit = 'conduitipid'
field_date = 'bv_cable_placed_planned'
# Populate a dictionary ...
Top 50 recent answers are included
