New answers tagged python
2
The data appears to be in "Nexrad Archive" format.
I installed PyART for Python 2 from git and its dependencies via pip:
https://arm-doe.github.io/pyart/
I could then run radar_info on one and get this:
$ radar_info KBBX20190819_012831_V06
## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar ...
0
In this case, you have to change the map, no the scale, you can use this:
map = layout.itemById("principal map")
map.setScale(12500)
if you don't have a itemByID, the other option is this:
layout = project.layoutManager().layoutByName("RAMMS_Karte_layout_de")
page = layout.pageCollection()
for i in items:
if isinstance(i, QgsLayoutItemMap):
...
0
For corroborating what kind of result your code produces, I created a qpt file (for this template https://raw.githubusercontent.com/qgis/QGIS/master/tests/testdata/template.qpt) and I loaded it in QGIS 3.8 with following code:
project = QgsProject.instance()
composition = QgsPrintLayout(project)
document = QDomDocument()
template_file = open('/home/zeito/...
1
Like BERA commented it is WKT.
In FME you could use a GeometryReplacer to (re)create the geometry of the feature based in the WKT string.
1
The error message is not very helpful unfortunately.
You are missing the required library LASzip on your system as well as the laszip.exe tool.
For Windows you can get the LASzip DLL file from the archive the LAStools project provides at https://lastools.github.io/download/LAStools.zip
You can find it in LAStools/LASzip/dll/ in said archive.
For quick ...
0
Finally cracked this (I am now using Win10, ArcPro 2.2.4 and R 3.6.1, but that shouldn't matter):
library(reticulate)
arcpy3_dir <- 'C:/Program Files/ArcGIS/Pro/bin/Python'
system2(file.path(arcpy3_dir, 'Scripts', 'proenv.bat'))
Then you can
import('arcpy')
and proceed, or use {python} chunks in an R markdown document e.g.
```{python}
import arcpy
/*...
0
The answer above from cm1 works well with PyCharm 2018. However, with PyCharm 2019 (2019.1 and 2019.2) it won't work anymore.
See also
PyCharm 2019 is not working with QGIS anymore
0
Assuming Address is unique to a customer (may be user customer ID instead) and that there can only ever be water or sewer as service types then a simple aggregate by customer counting on type will give you a 1 or a 2 for each customer address. You don't state in your question which GIS system you are using so assuming you are using ArcMap then this is the ...
0
Here is one way to do it (the simplified way ) :
let's suppose that in each shapefile you have a unique column called index, if you don't you can create it by reindexing the DataFrame
lines.reset_index(inplace = True)
grids.reset_index(inplace = True)
1 - first apply a spatial join between both of your shapefiles, with operation within to avoid extreme ...
3
The first warning means that the value of GeoASCIIParams tag is not read as it was written because the original image is having NULL character in the value of the tag. NULL can be used as a delimiter between strings http://freeimage.sourceforge.net/fnet/html/A633E9A9.htm but obviously GDAL takes just the first string.
The second error means that the writer ...
2
You can suppress the warnings (as long as you're sure there's no real issue with your data) with gdal.PushErrorHandler('CPLQuietErrorHandler'). If you do this any errors will also not get printed, so make sure you tell GDAL to raise a Python exception when an error occurs with gdal.UseExceptions().
E.g.
# Stop GDAL printing both warnings and errors to ...
0
I would recommend looking into the area of use for each of the CRS you are using in this example.
For the two options you have suggested for the input CRS, the bounds (min_lon, min_lat, max_lon, max_lat) are (113.76, 22.13, 114.51, 22.58).
>>> from pyproj import CRS
>>> CRS("ESRI:102140")
<Projected CRS: ESRI:102140>
Name: ...
0
You need to rescale the transform as well as the data, i.e.
from rasterio import Affine
from rasterio.enums import Resampling
scale = 10 # Reduce/upscale resolution scale factor
t = src.transform
# rescale the metadata
transform = Affine(t.a * scale, t.b, t.c, t.d, t.e * scale, t.f)
height = int(src.height / scale)
width = int(src.width / scale)
band = ...
0
I like to use gdal for these kind of operations. If you have access to the tools via the command line, it is very easy to resample your imagery at different resolutions like this example that would change "file1.tif" to a new file called "file1_0.5m.tif" which has been resampled at a resolution of 0.5m x 0.5m.
You can use gdalinfo to compare the resolution ...
1
I hope I understood your problem correctly (although still don't know why you added timestamps at the end if your csv clearly shows the time being two string columns).
First of all, you have to define an origin in space for your tensor. This origin consists of a lon, lat pair of coordinates that represent the upper left corner of your array. Let's assume ...
0
To add to what has been said by Michael, I would recommend computing the surface roughness of your DEM using the Rumple index or a simliar metric. You can also perform the roughness estimate on the point cloud itself as long as the ground points have been classified.
You may be able to classify the type of debris you are interested in based on the ...
2
One approach is to use the SQLite SQL dialect. It seems to work at least with a point shapefile where I digitized a few points, some of them in the same location.
ogrinfo -dialect sqlite -sql "select geometry, count(*) from duplicate_points group by geometry" duplicate_points.shp
Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 4
Extent: (271....
0
My idea is to order point coordinates by its longitude and latitude in this way:
Sort by longitude, in for loop iterate over the coordinates while longitude ascends;
Then iterate while latitude descends;
etc...
Then merge lists with coordinates and create polygon.
Not pretend this is the most elegant solution, maybe you send your data so I will test my ...
1
The simpliest way I think is to use boundary method in Shapely lib. When using
from shapely.geometry import Polygon, MultiPolygon
Polygon(feature['geometry']['coordinates'][0]).boundary
it returns
LINESTRING (30.916671 55.61667, 31.166671 56.91667, 32.550004 57.350003,...)
In case Multipolygon it obviously returns Multilinestring:
multy = MultiPolygon([...
1
The objects (Most of them, if not all) in GEE are immutable so you'll have to do
...
newList = newList.add(mosaic)
...
0
Ok, looks like I figured out why this was happening. In my post above I wrote that I was using version 3.7.3, but I had also tried running QGIS with version 3.6.8 using pyenv. I kept getting the same errors using 3.6.8 as well.
What solved my problem was to install the version of Python (3.6.8) that came with the QGIS (3.8) installer and the error messages ...
1
Try:
buffer_gdf = gpd.GeoDataFrame(geometry=buffer_polygs)
Another way using from_postgis copied from the answer to Geopandas PostGIS connection:
import geopandas as gpd
import psycopg2 # (if it is postgres/postgis)
con = psycopg2.connect(database="your database", user="user", password="password",
host="your host")
sql = "select geom, x,y,z from ...
1
You can convert x,y to col, row coordinates with gdal.InvGeoTransform, and gdal.ApplyGeoTransform. You can reproject coordinates with osr.CoordinateTransformation.TransformPoint.
E.g.
from osgeo import gdal, osr
ds = gdal.Open(someraster)
gt = ds.GetGeoTransform() # Geotransforms allow conversion of pixel to map coordinates
crs = ds.GetProjection()
...
3
No, it's not a bug. 5487115.2521567 is correct.
The GeoTransform doesn't give you the lower left, it gives you the origin. Which is the upper left. You can calculate the lower left from the origin + y pixel size * no. rows:
yllcorner = dem_transform [3] + dem_transform[5] * dem_file.RasterYSize
print(yllcorner)
5480455.0079735
0
As explained in comments, and demonstrated in the other answer (but not fully explained there), the solution was to make sure that each line is created with the same spatial reference as the polygons.
(I had actually tried this earlier, and got an error stating that spatial_reference is a read-only attribute. Turns out I had specified it incorrectly, and ...
0
This answer doesn't read the 2 shapefiles in using python but should do the job.
Load both shapefiles into QGIS
Run the "Line Intersections..." tool (Vector > Analysis Tools > Line Intersections)
Select one of your shapefiles as the 'Input layer' and the other as the 'Intersect layer'
Run the process to create a new point shapefile
Open the attribute table ...
2
Method "within" works fine, because this script
import arcpy
import itertools as itt
lines="LINES"
pgons="PGONS"
d=arcpy.Describe(pgons)
SR=d.spatialReference
g=arcpy.Geometry()
curT=arcpy.da.InsertCursor(lines,"Shape@")
with arcpy.da.SearchCursor(pgons,"Shape@") as cursor:
for row in cursor:
pgon=row[0]
gList=arcpy....
0
I kind of gave up on Minkowski's sum as it is too difficult for me to understand C/C++ code. As I was browsing one more time the excellent documentation of CGAL for 2D Minkowski Sum, I realised offsetting convex polygons or triangles should be easy enough.
So I resorted to triangulate my polygon using a readily available algorithm ('3d:tessellate' in ...
1
You can do this by combining min() and apply() functions of pandas
DF['associated'] = DF['vms_distance'].apply(lambda x: 'vms' if x == DF['vms_distance'].min() else 'no_vms')
0
Going to answer my question based on the reply found here. Apparently it is a bug, and the current workaround is to use:
myGeod.geometry_length(np.array(shapelyObject.coords))
instead of
myGeod.geometry_length(shapelyObject)
Will update when a final solution is available.
1
RQGIS only works with QGIS2, and I assume that you have already installed QGIS3. Hence, use RQGIS3, you can find it under https://github.com/jannes-m/RQGIS3. So far, it is not on CRAN due to an RStudio issue under Linux/Mac.
1
I think one of the most evident advantages is indexing. Consider the following example where you have data and both longitude and latitude stored in 2D arrays:
data = np.random.randint(100, 1000, size=(4, 4))
data
[[176, 479, 713, 973],
[992, 259, 969, 355],
[182, 139, 633, 938],
[761, 911, 124, 855]]
x = np.linspace(-76, -74.5, 4)
y = np.linspace(-5.0,...
0
I have encountered the same shift issue. This shift is not limited to GDAL/QGIS. It is also present in ArcMap, MapInfo and Global Mapper (maybe gdal is used in the backend...). A work around is to convert the GeoTiff to an ERS, fix the easting and northing in the .ers file (open with a text editor) and convert the ERS back to a GeoTiff. The conversion ...
1
Off course after I've asked the question I know the answer. After using Searchcursor I still have to iterate through it.
Changed the last line in the else to this and it works:
with arcpy.da.SearchCursor(inl, '*', where_clause=ftext) as feats:
featureIterator = []
for f in feats:
...
0
Here's an idea, creating a function which automatically creates a polygon with the area of 1 hectare for you according to the following steps:
Input the central coordinate
Reproject to a projection using metres as units (e.g. WebMercator)
Create a Square with the area of exactly 1 hectare
Project back to WGS 84
Here's an example implementation outputting a ...
0
I think your json is being interpreted as a string rather than as a json object.
Try importing the json library (part of the standard Python library) and then loading the string as a json object. For example:
import json
json_string = """
{
"pipeline":[
"D:/Lidar_collect_3/output/output_1.las",
{
"type":"filters.outlier",
"...
1
I believe I have now resolved this.
The polygon in the .kml file is missing a closing point (i.e. the same as the first point). Adding the first point at the end of the coordinates tag, importing to Google Earth Pro and measuring the area produces 955km^2, much closer to that produced by the python code.
2
A couple of things here:
Firstly don't set gdata to None until you are COMPLETELY done with that raster, data will return a broken dataset if you do. Please see for a more detailed explanation: https://trac.osgeo.org/gdal/wiki/PythonGotchas
Next thing uninstall gdal using pip then reinstall it using the .whl file from here:
https://www.lfd.uci.edu/~...
1
Sort by time and drop duplicate mmsis keeping last:
newdf = df.sort_values('time').drop_duplicates('mmsi', keep='last')
0
You can get the coordinates for each item and calculate how much of area of interest in covered by a particular item. I am using Shapely and Proj to do that in the following code block. Once I know which items I would like to download I use Planet's clip and ship API to only download the clipped areas.
from planet import api
import os
import json
import ...
2
The correct usage of RateLimiter for reverse geocoding would be:
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter
geolocator = Nominatim(user_agent="application")
reverse = RateLimiter(geolocator.reverse, min_delay_seconds=1)
location = reverse((50.6539239, -120.3385242), language='en', exactly_one=True)
print ...
0
If you want to show only two columns, you can make use of the setHorizontalHeaderItem() method for your model and then fill each row with values for both columns. The following is an example which could be run from the Python Console but you could adapt it for your plugin:
from PyQt4.QtGui import QStandardItem, QStandardItemModel, QTableView
layer = iface....
3
Something like this?
In the Shapely manual you can find published code on how you could generate images from geometry objects (that could be constructed by WKT).
linestring.py
from matplotlib import pyplot
from shapely.geometry import LineString
from figures import SIZE
COLOR = {
True: '#6699cc',
False: '#ffcc33'
}
def v_color(ob):
...
1
Here is a function to convert a QgsRasterLayer to a numpy array without GDAL through using the block method of the QgsRasterDataProvider(link):
from numpy import array
def convertRasterToNumpyArray(lyr): #Input: QgsRasterLayer
values=[]
provider= lyr.dataProvider()
block = provider.block(1,lyr.extent(),lyr.width(),lyr.height())
for i in ...
0
Convert all KMLs within a directory to SHP, zip by name, and move to another output folder? [closed]
to output to the output folder, add folder name like this
for %f in (*.kml) do ogr2ogr -f "ESRI shapefile" "output\%f.shp" %f
Shapefile is a multi-file format, and OGR updates the files in a non-sequential way. You have to zip it afterwards.
Merging can also be done using -append,-update. Refer to this
1
You are over complicating this a bit. Certain types of data are suitable for interpolation whereas other data is more well-suited for simple binning. In the lidar realm, elevation and intensity are commonly interpolated but, other lidar attributes (class, return number, time-stamp, strip-id) are simply binned to a pre-existing grid. From an analytic ...
0
Updating my question as I found a solution for my use-case.
I used folium for plotting. I also had the lat-long cords of the same UTM cords and used them. Matplotlib was not used.
import folium
from pandas import read_csv
df = read_csv('my/location/data')
lat = df['latitude'].values
lon = df['longitude'].values
center_lat = np.mean(lat)
center_lon = np....
0
Here's another way to do it:
import geopandas as gpd
import numpy as np
# load an example polygons geodataframe
gdf_polys = gpd.read_file(gpd.datasets.get_path('nybb'))
It looks like the following:
# find the bounds of your geodataframe
x_min, y_min, x_max, y_max = gdf_polys.total_bounds
# set sample size
n = 100
# generate random data within the bounds
...
0
In addition to installing GDAL, you also need to install the GDAL Python bindings. You do this with the command
pip install gdal
I don't know much about your MacOS environment and how you installed Python, so if this does not work, you would need to adjust your environment settings first to make pip work properly.
I am also not sure why you are trying to ...
1
After struggling one week I'm able to configure tilecache in Linux o.s. I also tested alongwith OL3 and working fine. The main problem is that the it is not reading tilecache.cfg files from correct path. So you have to give your absolute path for cfg files wherever required like in tilecache.cgi and services.py files. For mapserver service the mapscript ...
Top 50 recent answers are included
