texts: PyQGIS autoexport map canvas to PNG file with choosing location option

texts:
Following the query here:
QGIS capture screenshot of visible map (PyQGIS or other alternatives)
I would like to have the Export PDF option confined to clicking the button on my keyboard.
By using the solutions available here:
https://datatofish.com/screenshot-python/
I tried to do something like this with my code:
 import keyboard

 # if keyboard.is_pressed('Alt+S'):  (optionally, although I can still use the Run function in my Python console)

 file_path = filedialog.asksaveasfilename(defaultextension='.png')
 iface.mapCanvas().saveAsImage(file_path)

But I am getting errors i.e.

Traceback (most recent call last):   File "C:\PROGRA~1\QGIS3~1.12\apps\Python37\lib\code.py", line 90, in
runcode
exec(code, self.locals)   File "", line 1, in    File "", line 5
file_path = filedialog.asksaveasfilename(defaultextension='.png')
^

IndentationError: expected an indented block**
How can I save my screenshot with choosing and file name option?
I used another approach:
    iface.mapCanvas().saveAsImage(QFileDialog.getSaveFileName(mainWin, "Save
    file", "", "PNG Format (*.png)"))

but I am getting an error here:

Traceback (most recent call last):   File "C:\PROGRA~1\QGIS3~1.12\apps\Python37\lib\code.py", line 90, in
runcode
exec(code, self.locals)   File "", line 1, in    File "", line 1, in  TypeError:
QgsMapCanvas.saveAsImage(): argument 1 has unexpected type 'tuple'


texts: Treating latitude/longitude as a pair of polar coordinates

texts:
I have latitude longitude coordinates covering the whole globe that I wish to perform arithmetic operations (mean, interpolation) on without concern for ambiguous points and wrap around and without using geographic coordinate specific implementations.

In order to do this, I plan to convert these coordinates to a cartesian form.

i.e. for the position

29.7604° N, 95.3698° W

r=1 φ=29.764°, r=1 φ=95.3698°

I intend to then convert these polar coordinates to Cartesian coordinates

x=cos(29.764°) y=sin(29.764°), x=cos(95.3698°) y=sin(95.3698°)

In this Cartesian form can I then safely perform operations such as mean and linear interpolation?

I found this question/answer that seems to imply that the Cartesian form is insufficient, and that I would need a geographic specific representation: Averaging Cartesian ECEF coordinates then converting to WGS84 without latitude drift

But then I found this question/answer that seems to imply my strategy might just work: Computing an averaged latitude and longitude coordinates

What is the correct way to handle this?


texts: Read multiple layers of KML file using R

texts:
I have a KML that 3970 layers. I am trying to read it into R using the following script:

library(rgdal)
lyr <- ogrListLayers("file.KML")
mykml <- readOGR("file.KML","Layer1")

This works but what if I want to read all layers; something like:

mykml <- readOGR("file.KML",lyr)

Is there a way to accomplish this?


texts:
I could not open/read the following .img file.

The document below says it is a raster file.

http://edc2.usgs.gov/1KM/land_sea_mask.php

I tried to open it using many ways such as Q-Gis software, and Gdal translate, but could not.

Gdal translate reports error that unrecognized file..

That strange file is here:

http://edcftp.cr.usgs.gov/pub/data/avhrr/mask/wvsdn_mask_world.img.gz (3.54 MB).
I extracted it using winrar and the extracted file is 662 MB.I tried to open the extracted one.

Is someone familiar with this file?


texts: Counting points in shapefile within polygons of another shapefile using ArcGIS Desktop?

texts:
I'm creating a map with two shapefiles:

A polygon of counties
Points of cities

I used the function Intersect to group the cities in each county.

I'm trying to get the number of cities in each county to show up in the legend.

For example, in county A there are 16 cities, in county B there are 7 cities and so on.

I'm trying to create a choropleth map that shows the number of cities in each county.


texts:
When I look at the map "Prefeitura de Manaus (Gestão de Obras - SEMINF)" on ArcGIS online, it shows perfect alignment between the lot lines (in blue) and the satellite basemap:

However, when I look at the same map in the Explorer for ArcGIS app, there appears to be a distortion in the projection of the lots onto the satellite map, with all lots being misalligned by a few meters:

Does anyone know why the projections differ between the app and online and if there is a way to correct the projection in the app?


texts:
Is it possible, and if so how, to set the maximum scale a WMS layer can be viewed?

I know you can set it in a rule in the SLD style, but from the client you can always(?) change the style and thereby have other max scale limits. And that's not what I want.

I want to on the server side set that this layer should not be possible to view for scale larger than 1:1000.

Or is there any other good solution to accomplish this?


texts:
I have part of a LANDSAT image, this one: http://ge.tt/80uEl0n2

Using R, I would like to open this image, so I did:

B4<-raster(LE07_L1TP_219076_20021029_20170127_01_T1_B4.tif)

But it gave me this warning:

  Warning message:
  In .rasterFromGDAL(x, band = band, objecttype, ...) :

  This file has a rotation
  Support such files is limited and results of data processing might be wrong.
  Proceed with caution & consider using the "rectify" function

Then, using the rectify() function it gave me:

  Error in if (value[1] != nrow(x) | value[2] != ncol(x)) { :
  missing value where TRUE/FALSE needed
  In addition: Warning message:
  In dim<-(*tmp*, value = c(nr, nc)) :
   NAs introduced by coercion to integer range

How can I solve this to work with my tif image? At least, when I call B4, it gives me:

  class       : RasterLayer

  rotated     : TRUE

  dimensions  : 505, 872, 440360  (nrow, ncol, ncell)

  resolution  : 4.662376e-11, 1.623406e-10  (x, y)

  extent      : 259785, 285945, 7476655, 7491805  (xmin, xmax, ymin, ymax)

  coord. ref. : +proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

  data source : D:\L7Referencia\LE072190762002102901T1-
  SC20171002122402\LE07_L1TP_219076_20021029_20170127_01_T1_B4.tif

  names       : LE07_L1TP_219076_20021029_20170127_01_T1_B4

  values      : 0, 255  (min, max)

What I would like to do after this is open the other band and proceed calculating a vegetation index, like NDVI. But I can't even open properly them.


texts:
Trying to trasnform GCS_Bessel_1841 PCS_Trasnverse Mercator (central meridian 18) to WGS1984 using "Project tool" in ArcGIS 9.3. I tried to do it in several different ways, every time from the beginning.

In project dialog box under Output coordinate system I choose

a) Select predefined PCS. I choose WGS PDC. Nothing happens.
b) Select predefined GCS. I choose World/WGS_1984. Nothing happens.
c) Import from an existing dataset. Nothing happens.
d) I go to Layers, Transformations, select "Convert from GCS_Bessel_1841 to World/WGS_1984 (from the drop down menu). I select "using New Geographic Transformation" and choose position vector with seven variables which I got from here

http://georepository.com/transformation_3964/MGI-1901-to-WGS-84-4.html.

Nothing happens. The polygons remain the same, and the PCS is always the TM I started with.

Please help. Obvious ideas are welcome. As are a few short sentences behind the rationale of the system. I know what a PCS and GCS and datums are :) The country is Croatia.


texts: Fixing ArcScene Menubar that disappeared?

texts:
My ArcScene menubar is disappeared. I am trying multiple time to fix it but failed.

Can anyone help me to fix it?


texts: Split a layer using polygons from another layer

texts: I have a landuse layer that I would like to separate into different administrative units. What is the fastest way to achieve this in QGIS?


texts: Converting Pandas DataFrame to GeoDataFrame

texts:
This seems like a simple enough question, but I can't figure out how to convert a Pandas DataFrame to a GeoDataFrame for a spatial join?
Here is an example of what my data looks like using df.head():
    Date/Time           Lat       Lon       ID
0   4/1/2014 0:11:00    40.7690   -73.9549  140
1   4/1/2014 0:17:00    40.7267   -74.0345  NaN

In fact, this dataframe was created from a CSV so if it's easier to read the CSV in directly as a GeoDataFrame that's fine too.


texts: What to do with the intercept output from regression analsyis?

texts:
In ArcGIS I use global Ordinary Least Squares Regression (OLS) models and Geographically Weighted Regression (GWR) models to analyze the relationships between independent variables and my dependent variable.

What ArcGIS also includes in the results is an intercept term. I know basically the meaning of the intercept, but what I am unsure is - has the intercept to be interpreted and if yes, how exactly?

For GWR I get for each local unit an intercept - is it necessary/important/good to see the local variations of the intercept?


texts: Network Analysis: Solve Route - No Solution

texts:
I'm using ArcGIS 9.3.1

In my network there are isolated network(not connected to main network) and restricted(non-traversable network element position)

When I plot 1 point close to a isolated or restricted road. It will snap to it and return me no solution found or restricted network.

May I know is there a simple solution to ask the solver to choose another line/point if it can't find a solution?

In arcgis 9.3.1 there is no "Exclude restricted portion of the network" and I don't wish to use the query builder to ignore FType that carry restricted options.

So the main objective is to make the solver pick another line/point to snap on to.


texts:
Short explanation of my problem:

On the picture two polygons are shown. The highlighted (yellow) polygon is above the green polygon. Now, two operations are performed.

yellowFeature.geometry().overlaps(greenFeature.geometry()) --> False
greenFeature.geometry().contains(yellowFeature.geometry()) --> True

The results of the two functions were to be expected. Now the yellow geometry is divided into two halves.

Again, the two operations are performed.

yellowFeature.geometry().overlaps(greenFeature.geometry()) --> True
greenFeature.geometry().contains(yellowFeature.geometry()) --> False

The result I did not expect. I expected the same results as before the cut operation. This behavior results in another program to a major problem, since the geometry is correct and the yellow geometry does not overlap.

How can I change this behavior? Does it have anything to do with the precision (decimal)?

My settings: Topological editing is enabled. Snapping tolerance 5px.


texts: Connecting two coordinates on Google Earth Engine

texts:
I have about 500 flight data, where it has been separated as two csv's file (departure and arrival), in which each csv has 500 coordinate (lattitude and longitude).
The data looks like this

time_s
lat_deg
lon_deg

1
2.0
30
-90

2
4.0
31
-84

3
6.0
35
-89

4
8.0
35
97

5
10.0
35
49

6
12.0
35
24

7
14.0
35
78

8
16.0
35
49

...
...
...
...

500
1000.0
22
45

I want to plot this coordinates on GEE, where each coordinate is shown as point, and then the flight data (for example the first row(1) of departure and arrival is connected as a line.
So far I managed to show the point, as it is looked like this

Is that possible to connect each flight data (departure and arrival) as a line


texts: How to read the error when writing vector files in QGIS 3.0

texts:
Upgraded to QGIS 3.0 and now I cant read writer errors. The following (post 2to3) code doesnt print the error to the console:

CRS=QgsCoordinateReferenceSystem().fromEpsgId(4326)
error = QgsVectorFileWriter.writeAsVectorFormat(vl,fp, "system", CRS, "ESRI Shapefile")
if QgsVectorFileWriter.hasError(error):
    print("uh oh an error")
    print(QgsVectorFileWriter.errorMessage(error))

rather than reporting the error I get:

TypeError: QgsVectorFileWriter.hasError(): first argument of unbound method must have type 'QgsVectorFileWriter'

How can I modify this so that the write error is reported?


texts: Modis Reprojection tool, BAD or missing Total_BANDS values

texts:
I am trying to convert Modis terra Level 2 data 3KM swath i.e "MOD04_3K" from HDF to GeoTiff. I am using Modis Swath tool. When I input my file it gives an error saying "BAD or missing Total_BANDS values". When I input the file for Level 2 data 10 km swath it does not give any error.
Questions:

What is the problem with 3 KM swath data i.e giving "BAD or missing Total_BANDS values error"
When I try to Run tool for 10 Km swath it asks for Geolocation file. I have no idea about Geolocation file. It does not come with the downloaded MODIS file.

Or if u can suggest any other tool for reprojection of MODIS data and the conversion from HFD to GeoTiff


texts:
I would like to calculate a summarising column in QGIS. I have a couple of columns lets say x1, x2, ..x12 with string entries. Now I would like to count the entries over all columns i.e. x1 is not null then 1, x2 is not null then 1+1=2, x3 is not null 1+1+1=3,  x4 is null 1+1+1+0=3 etc. For some columns I would like to add more e.g. x5 is not null then 1+1+1+0+4= 7.

I tried with case when  "x1" is not null then 1 when "x2" is not null then +1 when "x3"  is not null then +4 when...

Is there an easy way to do?


texts: Edit field schema of web layer in ArcGIS Pro

texts:
I would like to edit the field schema of a web layer in ArcGIS Pro via Design --> Fields but I get the error message "the table or feature class corresponding to this view is read-only".

To clarify, I opened a web layer (Hosted Feature Layer) I have created previously in ArcGIS Pro (Version 2.5.0) using a portal connection. I checked the Editing Status and it says my layer is Editable. However, when I try adding new fields to the layer it says it's read-only. See image below.

Does anyone know what this means and how I can edit the field schema? Note, I have admin rights.


texts: How to find elevation highest and lowest for a state or city

texts:
I am trying to find out the highest and lowest elevation for these types of geographical areas:

continent (lots of ways to get info)
country (lots of ways to get info)
state (for the US states no problem, but what about other countries' federal states?)
city

Perhaps there a web service that gives such kind of information. Especially state and city I have problems finding elevation info on.

If not, what is my alternatives? Maybe I need to do it "manually"? I found via Google SRTM 90m Digital Elevation Data, is this a good way to get the elevation data I need? But then how could I use this to find it for a state or city?


texts: i have tried using the 'Multiple Layer Selection' plugin but it only offers the selection in triangle .. is there a way to have the same effect but with a polygon instead ?


texts: Line dash not working OL3 (Openlayers 3)

texts:
I want to display  a dashed line in a vector layer, as the it's shown in the Documentation there is a propriety called lineDash , but it didn't work for me here my code  :

 var iconStyle_line = new ol.style.Style({//  My style definition
          fill: new ol.style.Fill({
          color: 'rgba(221, 99, 0, 0.5)'
          }),
          stroke: new ol.style.Stroke({
          color: 'rgba(221, 99, 0, 0.5)',
          width: 3,
          lineDash: [4]
          })
     });

 var vectorSource_line = new ol.source.Vector();//  My vector source

 var vector_line = new ol.layer.Vector({

         source: vectorSource_line,
         style: iconStyle_line
     });


texts: v.surf.idw on grass return : no data points found

texts:
I'm a newbie with grass and I'm using grass (6.4.3) with the GUI... running Ubuntu 13.04 for what it's worth.

Here is my region :

projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      ** unknown (default: WGS84) **
ellipsoid:  wgs84
north:      871548
south:      871035
west:       176353
east:       176851
nsres:      1
ewres:      1
rows:       513
cols:       498
cells:      255474

Here are the five first lines of the files from which I created a vector map with v.in.ascii :

871207.125|176633.609375|322.54595947|1
871209.75|176633.84375|323.36795044|2
871207.4375|176623.8125|321.85391235|3
871213.75|176628.359375|322.46286011|4
871207.25|176658.4375|324.73287964|5

In table manager I have :

cat dbl_1 dbl_2 dbl_3 str_1
integer double precision double precision double precision character

I'd like to get a raster from this, so I'm trying to use v.surf.idw.

I cannot understand what to put in the column argument of this function.

When I put dbl_3, I get :
6423 records selected from table
WARNING : No record for point (cat = 6424)
WARNING : No record for point (cat = 6425)
0 points loaded
ERROR: No data points found

When I put double precision it is not working, I've tried every ways I can think of.

Thanks for reading me.


texts:
I just installed QGIS Wien 2.8.1, first the OSGEO4W Installer and later the standalone installer (32bits) and when I write import pyodbc on the Python console I get a message that the dll was not found, I did this on another computer (same bitness) and the same thing happens. However, when I use another computer which has QGIS 2.2 Valmiera (standalone installation 32 bit) and yet another which has 2.4 Chugiak (standalone installation 32 bit), import pyodbc works. I have not tried on 2.6.
All computers have Windows 7.


texts: How to make ArcCatalog recognize numbers in a csv/text file as numbers?

texts:
I am creating feature classes from GNS txt files in ArcCatalog (10.2). When I do this with the file for Bolivia (bl.txt), the feature class is created correctly. When I try to do this with the file for Argentina, the columns LONG and LAT are recognized as Text instead of Double, so I cannot assign them as X and Y fields.

The formats for both files are identical. Both ar.txt and bl.txt contain only numeric values in the columns LONG and LAT, with . as decimal separator. I cannot identify any differences between the two files to explain the different behavior. Is there a way to force ArcGIS to recognize these columns as double precision fields?

As a workaround I can import the txt  files in an Excel spreadsheet and create the feature class correctly from that, but I believe I shouldn't have to.


texts: Filtering (non-WFS) features in OpenLayers

texts:
Say I have a vector layer with a couple hundred features, and I just want to operate on a filtered subset of those features, such as those that have a State attribute of "OR".

Am I reading the OL docs correctly that the Filter classes are only for filtering WFS requests? Is there some other way of filtering ordinary features by a specific attribute (e.g. for zooming to the extent of just those features)?

Or do I have to just manually iterate over the features and populate an array of the desired ones? (Not a huge deal, by the way; just wondering if there's a more OpenLayers-y way of doing this).


texts: PostGIS - Find the "median" distance between two group of geometries depending on polygon area

texts:
I want to evaluate land distribution and concentration regarding to road network. Specificly, I'm trying to calculate the median distance between a shape of roads (multilinestring) and a shape of polygons (polygon) which will "split" the total polygon area in two. That is more or less a geographic median regarding a set of references geometries.

As I absolutly have no idea how to proceed, I tried some alternative approach.
I first calculate the nearest neighbor between each polygon centroid and road network and the minimal distance between those two. Then I calculate the median of those distances.

WITH nearest_distances AS
(SELECT
DISTINCT ON (polygons.id) polygons.id, roads.gid, ST_Distance(ST_PointOnSurface(geom), roads.geom) as distance
FROM polygons, roads
ORDER BY
  polygons.id, ST_Distance(ST_PointOnSurface(geom), roads.geom), roads.gid
)
SELECT
percentile_const(0.5) WITHIN GROUP (ORDER BY nearest_distances.distance)
FROM nearest_distances;

I a second time, I cross-check my median result with a query giving the area of polygons covered by a buffer on road network set on the result of the median.

SELECT
ST_Area(ST_Union(ST_Intersection(polygons.geom, ST_Buffer(road.geom, 220))))
FROM polygons, roads;

This is totally experimental as I needed quick estimation but I wonder if there any cleaner way to perform this. Idealy, I would split my polygons along a median line dividing polygon surface in two parts regarding the road network.


texts:
I'm working on a project for which I have GPS points and I would like to extract the sequence of roads related to those points. I know that this is called "map-matching" and that there is a load of literature about it.

Here is my question: is there a software / library that allows to easily convert a sequence of GPS points into a road sequence, given that I have the road network, so that I don't have to (re)implement the wheel?


texts: Adding new feature in QgsVectorLayer with openFeatureForm leads to segfault?

texts:
In a QGIS 2.18 plugin, I'm trying to create new features in an existing vector layer and I want to display the attribute form so the user can fill the attributes data.

The geometry is build before by parsing datas coming from websockets. Using layer.addFeature() works well for layers that does not have any fields constraints (eg: primary key).

As I need to write in Postgres/Postgis table with PK constraint, I want to use Feature Form to let the user fill attribute value and PK will handle by QGIS with default value. This way I also benefit from the settings that the user may have do for the layer fields (eg : combobox, default value, constraints, ...)

For that I'd like to use the openFeatureForm method of iface with my layer and my new feature like that:

feat = QgsFeature()
new_geom = QgsGeometry()# here I construct a geometry
feat.setGeometry(new_geom)
iface.openFeatureForm(self.layer, feat)

The attribute form show up as expected but when I hit the Ok button QGIS crash with segfault (core dumped) ...

Steps to reproduce in python console:

Select an editable vector layer
run this code snippet in QGIS python console

layer = iface.activeLayer()
 feat = QgsFeature()
 iface.openFeatureForm(layer, feat)
dialog show up!
click OK button --> Boom! QGIS crash with segfault

I guess something is wrong in passing feature to the form but i can't figure out how to make it works...

for the moment I'm using QGIS 2.18, it's planned to upgrade to QGIS 3.X but for now, I'm looking for a QGIS 2.X solution.


texts: QGIS calculate sum of values in more raster layers

texts:
I have cloud masks (rasters) of diferrent days (from sentinel 2a product) and I would like to calculate how many times was some area covered by clouds. My basic idea is put all layers on each other and calculate for each cell how many times value of cloud occured in this cell. Result should be a raster image with values of number of days under cloud for each cell.

How could I do that, what raster algebra should I use or should I convert rasters to vectors and made different processing ?


texts: Calculate Distance between 2 Shapely Points in Meters

texts:
I am trying to calculate the distance (in meters) between 2 lon/lng points of type shapely.geometry.point.Point using Python and Shapely. These 2 points are returned from PostGIS by SQLAlchemy/Geoalchemy2 as geoalchemy2.elements.WKBElement, which is then converted to shapely's shapely.geometry.point.Point using geoalchemy2.
However, when I try to use the distance function to calculate the distance in meters, the result is expected to be around 1490 m, but the code calculated it to be 0.013402616901217364.
What went wrong with our code? Thanks!
Note: Distance calculations will not exceed 30 km
from geoalchemy2.shape import to_shape

users = db_session().query(User).limit(2).all()
foo, bar = users[0], users[1]

print(type(foo.coordinates))               # <class 'geoalchemy2.elements.WKBElement'>
print(type(to_shape(foo.coordinates)))     # <class 'shapely.geometry.point.Point'>

distance = to_shape(foo.coordinates).distance(to_shape(bar.coordinates))
print(distance)                            # 0.013402616901217364


texts:
I am a game developer and I am not good in GIS, but it contains very useful data for our future project, my final goal is to generate height map for terrain render. I found a way to extract (convert) opened ASC file raster layer form 32 bit to 4 8 bits images, I changing min-max in layer properties, but I don't understand what units it uses

-24 to 74 ... I have no idea what that means.

P.S. I was searching forum for answer but did not find it, maybe I just don't know what I am looking for.


texts: SAGA dissolve tool keep only one column of data

texts:
I have installed QGIS 2.12 Lyon and I need to dissolve multiple polygons according to one column in the dbf file (my layer is a shapefile).

I have used SAGA "Polygon dissolve by attribute" algorithm and it works OK. The issue is that in my original file I have 7 columns, one of which is used to dissolve the polygons. The dbf file of the new layer has only 1 field - the one used to dissolve polygons.

How do I force SAGA to keep my other fields?

I just need one entry for each new polygon - information is repeated anyway.

Using the interface I could not see any option to select a method to retain information for other fields.


texts: Openlayers zooming and panning restrictions don't work

texts:
Consider the following example (I tried to keep it as simple as possible).

As the OSM map is in the Spherical Mercator Projection, the bounds (-1500000, 4000000, 5600000, 11000000) should cover most of Europe. I set these bounds both as the maxExtent and the restrictedExtent of the map.

Furthermore, I set maxScale to 60 M.

After launching this example, the map really centers to Europe and restricts panning to other parts of the World. But, it absolutely disregards the maxScale attribute and allows zooming far away.

Now, I try to remove the restrictedExtent. I should be able to pan the map freely. This is true. However, now the zoomToMaxExtent() doesn't work, although the maxExtent is still set!

How come neither maxScale nor zoomToMaxExtent() would work in this example? Am I just gravely misunderstanding the OpenLayers documentation?


texts: Append ArcGIS feature layer data using GDAL/OGR API

texts:
ArcGIS servers are configured with a maximum records count that they're allowed to return following a REST query. This maximum is usually 1000 records, which means that GDAL/OGR clients should implement an iterator to go through all the data in a layer.

My attempt to solve this problem was to use returnIdsOnly=true option to ArcGIS, which results in the list of object ids in the current layer. Next, I read these layers using a query on object id range, and this works fine as far as ArcGIS is concerned.

However, GDAL/OGR does not seem to have an easy way of concatenating results from different ArcGIS queries belonging to the same layer. Every time a data source is opened with OGRSFDriverRegistrar::Open function, a separate OGRDataSource is created with its own set of information.

I wonder if anyone could guide me on how to merge results of OGRSFDriverRegistrar::Open calls into the same OGRLayer object? I've tried using layer union operation, but it seems to be failing. Any hints?


texts: Getting a warning of "unexpected argument" for PyQGIS processing parameter in PyCharm

texts:
When specifying named argument for any type of QgsProcessingParameter, PyCharm complains of an Unexpected Argument warning. Although it doesn't cause any issue to the execution of the Processing algorithm in QGIS, it's kind of annoying when coding. I know I can turn off the inspection on this in PyCharm, but I hope there's some trick that I can use to make PyCharm understand the names (just as below: name, description and types). The complete definition of the class can be found here.

def initAlgorithm(self, config=None):

    self.addParameter(
        QgsProcessingParameterFeatureSource(
            name=self.INPUT,
            description=self.tr('Input layer'),
            types=[QgsProcessing.TypeVectorAnyGeometry]
        )
    )


texts: Qgis Dialog Not showing any Qwidgets Why?

texts:
i tried call another dialog form qgis plugin its show the dialog but not showing any Qwidget in that file
how to show them (if i call single ui file its showing fine )
Qgis Plugin Code:
from qgis.PyQt.QtCore import (QSettings, QTranslator, QCoreApplication,Qt)
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtWidgets import QAction,QDialog

from .ProgressDialog import PrograssDialog

def opendialog(self):
    Pdialog=PrograssDialog()
    Pdialog.ShowDialog()
def run(self):
    if self.first_start == True:
        self.first_start = True
        self.dlg = OpenNewDialog()
        self.dlg.button.clicked.connect(self.opendialog)
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:

PrograssDialog.py Code:
from qgis.PyQt import QtCore
from qgis.PyQt.QtWidgets import (QDialog,QPlainTextEdit,QScrollArea,QPushButton)
import sys
class PrograssDialog():
    def ShowDialog(self):
        try:
            self.Pd=QDialog()
            self.Pd.setWindowTitle("Script Execution... ")
            self.Pd.resize(500,500)
            self.ScrArea=QScrollArea(self.Pd)
            self.ScrArea.move(0,0)
            self.ScrArea.resize(500,300)
            self.TextArea=QPlainTextEdit(self.Pd)
            self.TextArea.move(0,0)
            self.TextArea.resize(500,300)
            self.Cancelbutton=QPushButton("Cancel",self.Pd)
            self.Cancelbutton.setGeometry(QtCore.QRect(360, 400, 93, 28))
            self.Pd.show()
        except Exception as msg:
            import sys
            tb = sys.exc_info()[2]

after execution plugin code it open dialog as shown in fig (No widgets) but i execute the PrograssDialog its showing all widgets.


texts: How to add geometry columns to a layer in QGIS 2.16?

texts: I want to add geometry columns to a shapefile. I use the tool called Export/Add Geometry Columns. The tool automatically creates a new layer with the input layer's attributes and the new geometry columns. In previous versions this tool (very usefully) allowed to add (or update) the columns to the existing shapefile instead of creating a new one. How can this be done in the 2.16 and 2.18 versions? I tried by selecting the curved arrow icon labelled ("Iterate over this layer") but the program crashes ("An error occurred while executing Python code").


texts: Sum area of features intersecting atlas feature in QGIS

texts:
I have AOI layer named "Clusters" which is set as the Atlas coverage layer for my project

I want to sum the area(Hectares) of specific intersecting features from my "Senwes_All_Fields_Alb_AND_IntStreek" layer in a text box in the print composer.

The "Senwes_All_Fields_Alb_AND_IntStreek" layer contains:

I would like to calculate the area of CATNAME_1 values intersecting the atlas feature using an expression

For example:

Horticulture = 500 Hectares
Non-pivot irrigated Annual Crop Cultivation / Planted Pastures = 2000 Hectares
Old Fields = 1000 Hectares
etc.

I have tried the aggregate function below, and got the TOTAL area of all CATNAME_1 values (Output preview: 347178.45651800046).

Now I would like to further break it down into the different values of CATNAME_1 as per example above.

Any ideas on how I might achieve that would be greatly appreciated.

aggregate(
layer:= 'Senwes_All_Fields_Alb_AND_IntStreek',
aggregate:='sum',
expression:="AREA_HA",
filter:=intersects( $geometry,geometry(@atlas_feature) )
)


texts: For periodic updates can a geodatabase table be replaced?

texts: I am pulling in CSV tables, that are updated frequently, into a file geodatabase.  Is there a way to overwrite the old table with the new one?


texts: Implementing genetic algorithm in ArcGIS Desktop?

texts:
I want to implement genetic algorithm to find optimum location for placement problem with ArcGIS.

Can any body explain it (step by step ) or any resources about that?


texts: Google Engine - Make a greenest pixel composite for Sentinel 2

texts:
Does anyone know how to do a composite of greenest pixel for Sentinel 2?

That is calculate the maximum NDVI for a Image collection, and then do a pixel-wise  image stacking of the original bands (that correspond to the image of the maximum NDVI).
There is already an algorithm for Landsat 8 named  .qualityMosaic()  , but it needs the quality band that is only available for Landsat.


texts: Combining a MapTool and a QGIS plugin toolbar?

texts:
This is a follow-up question to my question about how to build a QGIS plugin and link a textbox to the plugin toolbar (Adding textbox to QGIS plugin toolbar?).

In addition to that, I was wondering how to add functions of a Map Tool in order to start a different "mode". By mode, I mean something like a editor mode or the identify tool, where the meaning of mouse clicks, scroll event etc change to fit the current task at hand.

I found an answer to this GIS.SE question (Programatically check for mouse click in PyQGIS?), and tested its functions. It seems to be a good start and can handle a mouse click event, when activated with

tool = PointTool(qgis.iface.mapCanvas())
qgis.iface.mapCanvas().setMapTool(tool)

My questions are:

Both the PointTool(QgsMapTool) and plugin created with QGIS's standard Plugin Builder (and pb_tool) can handle various events. Which events should be handled by the MapTool vs by the Plugin? Is there a convention on this?

Also, how to link the MapTool to the Plugin (or embed the MapTool in the Plugin)? Is there a conventional way?

For example, if I click any button in the plugin toolbar, the mouse icon should change to a different shape. Currently, the mouse is still a magnifying class after I activate the MapTool.

How to properly register the tool when the plugin is activated and restore to the previous tool/mode when, e.g., the plugin toolbar is closed/hidden?


texts: How to listen group layer added layer in ArcObject 10.1?

texts:
I tried using the following code to listen to the layers loaded into MapControl. However, it is found that the next event is not triggered when the layer is added to the layer group. How do you listen to the layer group to add and delete layers?

mEvent = pAxMapControl.ActiveView.FocusMap as IActiveViewEvents_Event;
mEvent.ItemAdded += new IActiveViewEvents_ItemAddedEventHandler(Event_ItemAdded);
mEvent.ItemDeleted += new IActiveViewEvents_ItemDeletedEventHandler(Event_ItemDeleted);

private void Event_ItemAdded(object Item)
{
    if (Item is IFeatureLayer)
    {
        FeatureLayerAdded(Item as IFeatureLayer);
    }
    else if (Item is IRasterLayer)
    {
         RasterLayerAdded(Item as IRasterLayer);
    }
    else if (Item is ILasDatasetLayer)
    {
         LasLayerAdded(Item as ILasDatasetLayer);
    }
}


texts:
I work in Arcmap 10.3.1.

I have raster dataset - spatial reference 3413

My problem is that I can not send this data on a map.

Sequence:

Add raster dataset
Add basemap

How do I display the data on the card in the correct projection?


texts: arcpy taking longer and longer in a loop

texts:
I am looping through an amount of point arrays, creating data on the fly.

The features/arrays, are largely similar in size.  One thing I have noticed, is that when I start the process, the time it takes to perform, rises from about 8 seconds, to around 1 minute.  I can't see why this is happening, but you cna see the times slowly increase if you put a timer around the loop.

Is there a way to increase this performance? Is there a flush method?  I believe I am properly handling my objects and deleting them after use, but I cannot for the life of me see why it is happening?

Any ideas?

The code is refactored, so there may be some typos, etc.

def createGeom(geom,header,band,scratchDB,buffer, dist):

filetime = (str(time.time())).split(".")
outfile = "fc" + filetime[0]+filetime[1]
outpath =  scratchDB + "tmp.gdb/Polygon/"
print outpath
sCon = ["1", "16", "2.0"]
ts = header.split(" ")
ti = ts[6]
Current = ts[9]
LL = ts[2][2:5]
LU = ts[2][8:11]

outFeatureAggClass = outpath+outfile+"_agg"
outFeatureClass = outpath+outfile

arcpy.AggregatePoints_cartography(geom, outFeatureAggClass, dist)
arcpy.Buffer_analysis(outFeatureAggClass, outFeatureClass, buffer)

arcpy.AddField_management(outFeatureClass, "Name", "Text")
arcpy.AddField_management(outFeatureClass, "ID", "SHORT")
arcpy.AddField_management(outFeatureClass, "LL", "SHORT")
arcpy.AddField_management(outFeatureClass, "LU", "SHORT")
arcpy.AddField_management(outFeatureClass, "Tim", "SHORT")
arcpy.AddField_management(outFeatureClass, "Band", "SHORT")
arcpy.AddField_management(outFeatureClass, "Lvl", "TEXT")
arcpy.AddField_management(outFeatureClass, "Cur", "SHORT")

rows = arcpy.UpdateCursor(outFeatureClass)
for row in rows:
    row.name = outfile
    row.ID = 1000
    row.LL = int(LL)
    row.LU = int(LU)
    row.Tim = int(ti)
    row.Band = int(band)
    row.Lvl = sCon[int(band)]
    row.Cur = Current
    rows.updateRow(row)

del row
del rows

UPDATE

I am now running this in ArcMap, and it's taking significantly less time, to a factor or about 3 or 4; what took 30 mins, now takes 7.

What's that all about?

EDIT

ESRI have agreed to look at the issue. I'll park this for now, ready to update it with any progress.


texts: arcpy.mapping.RemoveLayer function not removing layer from mxd

texts:
I have a Python script that ends with a function to remove a feature layer, "Address Point", that is created at the start of the script. When I run the script it completes without any errors but does not remove the layer that is defined in the RemoveLayer function. If I list out the layers in the mxd after the script completes the layer I am trying to remove is listed with the layer name I am calling. Any ideas on why this function is falling to remove this layer?

    import arcpy
    import os

    arcpy.env.overwriteOutput = True

    # Make Layer
    arcpy.MakeFeatureLayer_management(r'Database
    Connections\10.103.2.145_BCDES_sa_ArcGIS_10.2.2.sde\BCDES.DBO.AddressPoint','Address Point')

    # Delete Feature class from file Geodatabase
    arcpy.Delete_management('C:\GIS\BlairMailer\Mailer.mdb\Address_Point')
    print 'Feature Class Deleted'

    # Select by attribute
    arcpy.SelectLayerByAttribute_management("Address Point", "NEW_SELECTION", "[PRINT_] = 'Y'" or "[PRINT_] = 'y'")
    print 'Features Selected'

    # Export selected feature to Geodatabase
    arcpy.FeatureClassToGeodatabase_conversion("Address Point",'C:\GIS\BlairMailer\Mailer.mdb')
    print 'Feature Class Exported'

    # Open Mailer
    os.startfile(r'C:\GIS\BlairMailer\Mailer.docx')

    # Calculates Null valve for selected features
    arcpy.CalculateField_management("Address
    Point","PRINT_","None","PYTHON_9.3")

    #Remove feature layer from map
    mxd = arcpy.mapping.MapDocument(r'C:\Users\achapman\Desktop\Addressing.mxd')
    for df in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
        if lyr.name.lower() == "Address Point":
            arcpy.mapping.RemoveLayer(df, lyr)

Revised my script and this is working now.

import arcpy
import os

arcpy.env.overwriteOutput = True
currentMXD = arcpy.mapping.MapDocument("current")
activeDataFrame = currentMXD.activeDataFrame

# Make Layer
arcpy.MakeFeatureLayer_management(r'Database Connections\10.103.2.145_BCDES_sa_ArcGIS_10.2.2.sde\BCDES.DBO.AddressPoint','Address Point')

# Delete Feature class from file Geodatabase
arcpy.Delete_management('C:\GIS\BlairMailer\Mailer.mdb\Address_Point')
print 'Feature Class Deleted'

# Select by attribute
arcpy.SelectLayerByAttribute_management("Address Point", "NEW_SELECTION", "[PRINT_] = 'Y'" or "[PRINT_] = 'y'")
print 'Features Selected'

# Export selected feature to Geodatabase
arcpy.FeatureClassToGeodatabase_conversion("Address Point",'C:\GIS\BlairMailer\Mailer.mdb')
print 'Feature Class Exported'

# Calculates Null valve for selected features
arcpy.CalculateField_management("Address Point","PRINT_","None","PYTHON_9.3")

# Open Mailer
os.startfile(r'C:\GIS\BlairMailer\Mailer.docx')

# Remove Address Point layer
layerList = arcpy.mapping.ListLayers(currentMXD, "", activeDataFrame)
for layer in layerList:
    if layer.name == "Address Point":
         arcpy.mapping.RemoveLayer(activeDataFrame, layer)

arcpy.RefreshTOC()


texts: ajax don't retrieve fully php request response

texts:
I'm facing challenge with ajax response in my code:

//  ********* PHP code (LoadData.php) ****************
<?php
include 'Connexion.php';

$requete="SELECT id, pkuid, province, nom, prenom, tel, nom_mammifere, qte_male, qte_femell, ST_AsGeoJSON(geom, 5) as geom FROM mammiferes ORDER BY nom";

$Connection = Connexion::getInstance()->getDB();
$res = $Connection->prepare($requete);
$res->execute();

$features=[];
while($row = $res->fetch(PDO::FETCH_ASSOC)){
    $feature =["type"=> "Feature"];
    $feature["geometry"] = json_decode($row["geom"]);

    unset($row["geom"]);

    $feature["properties"]=$row;
    $features[]=$feature;
    //array_push($features, $feature);
 }

$featureCollection = ["type"=>"FeatureCollection", "features"=>$features];
echo json_encode($featureCollection);
?>

This works fine and retrieve data using such url: localhost/myproject/loadData.php
Now I want to send the output of this file in my javascript file using ajax. Below is my ajax call:

$.ajax({
url:'loadData.php',
type:'POST',
success:function(response){
    alert(response);
},
error:function(xhr, status, error){
        console.log(xhr);
        console.log(status);
        console.log(error);
    }
 });

this retrieve just a part of the php result (for example 75 of 3500 rows) of this result is sometime not a valid json (jsonlint.com).
I don't know what is going wrong. Merci


texts: How to "Select By Location" in ArcGIS Server JavaScript API

texts:
Is it possible to perform the equivalent of "select by location" in the ArcGIS Server JavaScript API v4.x? The use-case is "find features in Layer 1 which are within x kilometres of features in Layer 2".
Here are the options I've considered:

a QueryTask with Query.geometry - this will query an entire layer, but only has the option to query against a single geometry object, not another layer
GeometryEngine.intersect - only works between 2 geometry objects, not 2 layers
The FindExistingLocations spatial analysis service, which seems to require special licensing which I don't have (I don't have access to ArcGIS Portal on this project)
Use the GeometryEngine.intersect approach in a JavaScript loop - which would require comparing every feature in both layers against each other (?)

I'm a little surprised that something so commonplace in a desktop GIS doesn't seem to be easily achievable in a web map. Am I missing something obvious?


texts: Comparing values across attributes of one feature in QGIS

texts:
I am using QGIS 2.18

I have a shapefile that contains data of each state of India, and the number of people per state following a particular religion.
 Each feature has the following attributes:
name, total population, total hindu, total muslim, total christian, (there are 8 fields on religion)etc etc.Taken from census data.

I want to do stylize my map in several ways that involve finding out the religion with the maximum people in each state.

So basically, the styling function I'd need performs something like this:

Case
when ~Hindu is max~ then color: color1
when ~christian is max~ then color: color2
....

I've been trying to see if I can compare across fields and then return the field name of which the maximum value belongs to. But I cannot find any function that does that.

Any tips on how to do this?

PS I have not studied POSTGIS,grass etc, so it would be great if I could be given answers that do not require these. I will be studying it soon though, so if it is inevitable to use that I will. But I want to see what I can do now.


texts: ArcObjects "Tool" add-in - Polygon not persisting in display while moving

texts:
I am creating a "tool" add-in for ArcGIS 10 that basically draws a polygon where user left-clicks on map display and, as long as mouse is left-clicked, drags it around. When left button is released, polygon is erased.

It should not be so complicated, but I'm struggling in making the polygon persist on screen while dragging. When mouse is clicked, polygon is displayed ok, but as soon as I start moving around it vanishes - but if I move the mouse fastly, I can see the polygon there (like a "flickering" effect). That is: polygon is drawn, but when mouse stops moving, it is cleared from screen.

If I remove the mxDoc.getActiveView().refresh() call, every polygon drawn during movement persist on screen. So may be something related to the correct refresh method to use? What am I doing wrong?

I'm just beginning my development with ArcObjects and Engine, so I might be missing something here.

Here is a piece of my code (Java):

@Override
public void mouseMoved(MouseEvent me) {
  try
  {
      if(me.getButton() != 0)
      {
          mxDoc.getActiveView().refresh();
          gabPolygon = new Polygon();
          gabPolygon = getPolygon(me); // my code - creates a Polygon object based on
                                       // mouse coordinates
          screenDisplay.startDrawing(screenDisplay.getHDC(), (short) esriScreenCache.esriNoScreenCache);
          screenDisplay.setSymbol(getSymbol()); // my code - sets polygon symbology
          screenDisplay.drawPolygon(gabPolygon);
          screenDisplay.finishDrawing();
      }
  }
  catch(Exception e){
      e.printStackTrace();
  }
}


texts: python add fields that may have NULL values

texts:
I would like to use python to calculate a field (short integer) using other fields (short integer) present in the table. In some cases these fields may have NULL values for one or more of the fields. For example, the new field name is 'total' and total is based upon adding field01 + field02 + field03. Some records (rows) have values for field01, field02, and field03, but some rows only have values for field1 and field2, then field3 is NULL.

I have some code that I believe will work, but I must not be implementing it correctly.

#calculate new field 'total' in feature class 'test'
inFeatures='test'
fieldName01='total'
exSummed = "stack(!field01!+!field02!+!field03!)"
codeBlock = "def stack(*args):
    return sum(filter(None, args))"

arcpy.CalculateField_management(inFeatures, fieldName01, exSummed, "PYTHON", codeBlock)

Here is the error I receive...

Thanks
-al

UPDATE: @Jason was correct, the triple double quote was my problem. I also modified my code based upon a suggestion posted here.

This is my updated code to reflect my solution:

import sys, string, os, arcpy
from arcpy import env

# Set workspace
env.workspace = "Q:/LandTrendr_04/lt_ancillary/landtrendr_file_geodatabase/templates_testing/LandTrendr_PIRO.gdb"

#=========================
#calculate new field 'total' in feature class 'test'
inFeatures='test'
fieldName01='total'
exSummed = "stack(!field01!,!field02!,!field03!)"
codeBlock = """def stack(field01,field02,field03):
    itemList = [field01,field02,field03]
    myList = [item for item in itemList if (item !=None)]
    return sum(myList)"""

arcpy.CalculateField_management(inFeatures, fieldName01, exSummed, "PYTHON", codeBlock)


texts: How to fix QGIS Plugins crashing with "NameError: name 'pyplugin_installer' is not defined"?

texts:
I recently switched to Windows 10 and therefor did a clean Windows and QGIS 2.10.1 Install. QGIS worked fine, I reinstalled all Plugins, that I used to have. No Problems.

Today I switched on the PC and wanted to open QGIS and was prompted with the following Error Messages. What bugs me is, that I had similar problems in my Windows 8 and QGIS 2.8 system.

The whole Python Plugin System seems to "crash". QGIS opens but all Plugins are disabled and I can't even access the Plugin-Module (where I can pick an install/uninstall plugins).

Could somebody help me with that please? A reinstall of QGIS 2.8 didn't help when I had a similar problem months ago.

The error messages (german installation):

  Fehler bei der Ausführung folgenden Codes: import pyplugin_installer

  Traceback (most recent call last):   File "", line 1, in    File
  "C:/PROGRA~1/QGISPI~1/apps/qgis/./python\qgis\utils.py", line 200, in
  loadPlugin
      msg = msgTemplate % (packageName, "', '".join(sys.path)) UnicodeDecodeError: 'ascii' codec can't decode byte 0xf6 in position
  18: ordinal not in range(128)

  Python-Version:
  2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

  QGIS-Version:
  2.10.1-Pisa 'Pisa', d20c5b7

  Fehler bei der Ausführung folgenden Codes:
  pyplugin_installer.initPluginInstaller()

  Traceback (most recent call last):   File "", line 1, in  NameError:
  name 'pyplugin_installer' is not defined

  Python-Version:
  2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

  QGIS-Version:
  2.10.1-Pisa 'Pisa', d20c5b7

  Python-Pfad: ['C:/PROGRA~1/QGISPI~1/apps/qgis/./python',
  u'C:/Users/Philipp B\xf6hme/.qgis2/python', u'C:/Users/Philipp
  B\xf6hme/.qgis2/python/plugins',
  'C:/PROGRA~1/QGISPI~1/apps/qgis/./python/plugins',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\bin\python27.zip',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\DLLs',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\plat-win',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\lib-tk',
  'C:\PROGRA~1\QGISPI~1\bin',
  'C:\PROGRA~1\QGISPI~1\apps\Python27',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\PIL',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\win32',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\win32\lib',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\Pythonwin',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg',
  'C:\PROGRA~1\QGISPI~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg',
  u'C:/Users/Philipp B\xf6hme/.qgis2//python', 'C:\Users\Philipp
  B\xf6hme\.qgis2\python\plugins\DigitizingTools\tools',
  'C:\Users\Philipp B\xf6hme\.qgis2\python\plugins\mmqgis/forms


texts: How to find road distance between two points in PyQGIS

texts:
I am working where I need to make python scripts to find distance between set of points using any predefined library function.

While searching for same I came across Road Graph plugin of C++.

Is there any similar library for Python in PyQGIS where I can find shortest distance among two points when road line layer is given.


texts:
This is a follow up question to this question.

I have a river network (multiline) and some drainage polygons (see picture below). My goal is to select only the headwater polygons (green).

With John's solution I can easily extract the river start points (stars).
However, I can have situations (red polygon) where I have startpoints in a polygon, but the polygon is not a headwater polygon, because it is flown though by the river. I only want the headwater polygons.

I tried to select them by counting the number of intersection between polygons and rivers (rationale: a headwater polygon should have only 1 intersection with the river)

SELECT
    polyg.*
FROM
    polyg, start_points, stream
WHERE
    st_contains(polyg.geom, start_points.geom)
    AND ST_Npoints(ST_Intersection(poly.geom, stream.geom)) = 1

, where poylg are the poylgons, start_points from johns answer and stream is my river network.

However, this take forever and i did not run it:

"Nested Loop  (cost=0.00..20547115.26 rows=641247 width=3075)"
"  Join Filter: _st_contains(ezg.geom, start_points.geom)"
"  ->  Nested Loop  (cost=0.00..20264906.12 rows=327276 width=3075)"
"        Join Filter: (st_npoints(st_intersection(ezg.geom, rivers.geom)) = 1)"
"        ->  Seq Scan on ezg_2500km2_31467 ezg  (cost=0.00..2161.52 rows=1648 width=3075)"
"              Filter: ((st_area(geom) / 1000000::double precision) < 100::double precision)"
"        ->  Materialize  (cost=0.00..6364.77 rows=39718 width=318)"
"              ->  Seq Scan on stream_typ rivers  (cost=0.00..4498.18 rows=39718 width=318)"
"  ->  Index Scan using idx_river_starts on river_starts start_points  (cost=0.00..0.60 rows=1 width=32)"
"        Index Cond: (ezg.geom && geom)"

So my question is: How can I efficiently query headwater polygons?

Update:
I added some sample data to my dropbox. Data is from south-west Germany. It's two shape files - one with streams and one with polygons.


texts: Creating a layer after using SelectLayerByAttribute_management

texts:
Is it possible to use an existent layer to do a selection by attributes without having to use MakeFeatureLayer_management first? Because I already have a layer with the same name as the feature class "Rede_Cbr_MM_2008_proj".

Only this code works for me:

import arcpy
from arcpy import env

# Set workspace
env.workspace = "S:/workspace/MSc_git/ArcGis/estagio_db.gdb"
workspace = "S:/workspace/MSc_git/ArcGis/Layers"
in_feature= "Rede_Cbr_MM_2008_proj"
tmp_layer="tmp_layer"

# # Set overwrite option
# arcpy.env.overwriteOutput = True

#MakeFeatureLayer variables
where_clause='TSYSSET NOT LIKE \'%Metro%\''

try:
    #Make a layer from the feature class
    arcpy.MakeFeatureLayer_management(in_feature,tmp_layer)
    arcpy.SelectLayerByAttribute_management(tmp_layer,"NEW_SELECTION",where_clause)

    # Write the selected features to a new featureclass
    arcpy.CopyFeatures_management(tmp_layer, "Rede_Cbr_MM_2008_WithoutMetro")
except:
    print arcpy.GetMessages()

But if I try to use my layer using only SelectLayerByAttribute_management I get this error:

Executing: SelectLayerByAttribute S:/workspace/MSc_git/ArcGis/Layers/Rede_Cbr_MM_2008_proj.lyr NEW_SELECTION "TSYSSET NOT LIKE '%Metro%'"
Start Time: Sat Jul 20 21:49:55 2013
Failed to execute. Parameters are not valid.
ERROR 000825: The value is not a layer or table view
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (SelectLayerByAttribute).

Isn't the SelectLayerByAttribute supposed to take a layer as 1st parameter? I don't know I do I have to do this kind of "layer on the fly" (tmp_layer) and then select.


texts: QGIS - 3D data set with stacked points at different elevations, is it possible to create data points at each elevation

texts:
I have a 3D dataset consisting of xyz and property. The data points are along vertical axis (boreholes), and at different XY locations (at other borehole location).  The data points are not all at the same elevation between each borehole location due to different soil intervals, and the property value is only provided for the Top of the interval.

My goal is to populate the borehole with properties between different soil Tops of interval, so that I can map properties at different horizontal elevation slices.

I found how to do this with Cross Join Query in MS Access but this creates very large data sets which can make the software crash as it provides every permutation possible.

Once that is done, the data needs to be filtered. So I am wondering if in QGIS, if there is a way to do something similar without having to generate this huge dataset by a Cross Join query.  For example it would resemble something like this in plain english: "show points at elevation 110 masl at each Borehole location, if no point at elevation 110, then pick next one up."


texts: How to change Leaflet Map panes layering order (z-index)?

texts:
Leaflet maintains that the Map panes elements contain all layers added to the map. The difference between many of the Map panes is strictly the z-index order of layering.

I would like to use a combination of lvector.CartoDB layers, which are essentially overlayPane layers, with TileLayer, such as GeoIQ Acetate-bg and labels.

This is the ordering of the elements as they are added the map:

tileLayer1 = new L.TileLayer();
map.add(tileLayer1);  // add first layer to map

cartoDBLayer1 = new lvector.CartoDB();
cartoDBLayer.setMap(map); // add second layer to map

tileLayer2 = new L.TileLayer();
map.add(tileLayer2);  // add third layer to map

What returns is a map with layers in this order:

tileLayer1,tileLayer2,cartoDBLayer1

tileLayer1 and tileLayer2 are situated in the HTMLElement: TilePane and cartoDBLayer1 is in HTMLElement: overlayPane.

Is there any way to force cartoDBLayer1 to render in the TilePane, such that it falls in order of the z-index that it is added to the map in...

i.e.

z-index[0]:tileLayer1
z-index[1]:cartoDBLayer1
z-index[2]:tileLayer2


texts:
I have a recurring problem in Civil 3D 2011.  From time-to-time, I have to attach to an image file (ie - survey plan) and draft on top of it.  In the case of the current work I'm doing, I am making revisions to an existing drawing set.  That said, we need to keep the original work, so I've saved the drawings to a new folder with a new name (indicating the date of my revisions).  I have set "relative pathways" on the images I've inserted, and the drawing (technically) should go fishing for the images starting at the new save location.  When I open the drawings however, the images are missing, and all I see is a link indicating where the images "should be".

In ArcMap, this type of problem is easy to fix.  Does anyone know how to fix this?

I know this is a GIS community, but I figure some of you must use CAD too.


texts: set bounds for pointLayer based on another layer size

texts:
I am creating a simple method for drawing Points on map, and have problem with full extent. I create a new point layer and add it to mapContent, and the problem is that it set bounds somewhere around 0, is it possible to set that bounds to be same as for another layer opened.
If I first open another layer I get this in console

Jun 22, 2016 6:26:54 PM org.geotools.map.FeatureLayer getBounds
WARNING: Bounds crs not defined; assuming bounds from schema are correct for CollectionFeatureSource:org.geotools.feature.DefaultFeatureCollection@48b99eac
Jun 22, 2016 6:26:55 PM org.geotools.map.FeatureLayer getBounds
WARNING: Bounds crs not defined; assuming bounds from schema are correct for CollectionFeatureSource:org.geotools.feature.DefaultFeatureCollection@48b99eac

and I get that every time I create a new point, just more for every new point

Here is code I use for this

    JButton create = new JButton("create");
    create.addActionListener(new ActionListener() {
        public void actionPerformed(ActionEvent e) {
            mapPane.setCursorTool(
                    new CursorTool() {
                        @Override
                        public void onMouseClicked(MapMouseEvent ev) {
                            createPointLayer();
                            createFeatures(ev);
                            for(Layer layer : map.layers()){
                                if(pointLayer == layer){
                                    i++;
                                    if(i>1){
                                        int index = map.layers().indexOf(layer);
                                        map.layers().remove(index);
                                    }
                                }
                            }
                        }
                    });
        }

    });
    toolBar.add(create);
    }
void createFeatures(MapMouseEvent ev) {
    double[] xy = ev.getWorldPos().getCoordinate();
    Point point = geometryFactory.createPoint(new Coordinate(xy[0], xy[1]));
    pointCollection.add(SimpleFeatureBuilder.build(pointType, new Object[]{point}, null));

    //create map layer event
    MapLayerEvent mple = new MapLayerEvent(pointLayer, MapLayerEvent.DATA_CHANGED);
    //create maplayer list event
    MapLayerListEvent mplle = new MapLayerListEvent(map, pointLayer, map.layers().indexOf(pointLayer), mple);

    mapPane.layerChanged(mplle);
    System.out.println(MessageFormat.format("Created Point: {0}", point));
}
/**
 *
 */
private void createPointLayer() {
    if (pointType == null) {
        pointFeatureTypeBuilder.setName("Point");
        pointFeatureTypeBuilder.setCRS(mapPane.getMapContent().getCoordinateReferenceSystem());
        pointFeatureTypeBuilder.add("the_geom", Point.class);
        pointType = pointFeatureTypeBuilder.buildFeatureType();
        pointCollection = new DefaultFeatureCollection(null, pointType);
    }

    pointLayer = new FeatureLayer(pointCollection, pointStyle);
    map.addLayer(pointLayer);
}

Do you have any suggestion, or to leave it like that and just exclude that layer when setting extent?


texts: How to mosaic two tif data by gdal_merge

texts:
my method is as follows(Some paths have been omitted.)

 enter code here gdal_merge.py F:\test_data\LC81210402013278LGN00_B1.TIF F:\test_data\LC81220402013285LGN00_B1.TIF -o F:\mo529.tif

but failed,The error is as follows：

How should I do?


texts: CSV to ArcGIS table, insufficient text field length

texts: I currently use arcpy.TableToTable_conversion in order to convert a CSV file into an ArcGIS file geodatabase table. This works when string length is less than or equal to 255 characters. Whenever 255 characters is exceeded the conversion fails. Is it possible to force the textfield in the table to be larger than 255 characters?


texts: Pricing for ArcGIS for Server 10.3 and Portal for ArcGIS?

texts:
I just had a meeting with an Esri distributor and learned that Portal for ArcGIS will not be an extension but instead will be a part of ArcGIS for Server beginning from 10.3.

We currently do not have either Server or Portal products and will invest on both of them. In current version (10.2) if we buy ArcGIS Server and the Portal Extension what we will pay for Portal will be as much as what we pay to Server.

Below is the question they cannot answer; Is there a possibility to have the price of the Server 10.3 product (which will include Portal as well) remain the same with 10.2?

Because if that will be the case, investing for both products now will be a total loss of money. But I do not really think they will double the Server price in 10.3 version since it will then be a problem for the customers that require Server but not Portal.

Anyone had similar past experiences with Esri and give me an idea how to proceed?


texts: QGIS Multi distance buffer units do not align with the main layer

texts:
I have a buildings shapefile of Beijing. I know the projection WGS_1984_UTM_Zone_50N
WKID: 32650 Authority: EPSG
So I would like to create multi distance buffer from a central point but when I do so, buffer of 1 m covers 3x size of the city. I change the projection of the layer when it was uploaded and checked all the settings to match. Tried to reopen the layer but still have the same problem. What am I doing wrong?


texts: OpenLayers view.getResolution() syntax error

texts:
I'm trying OpenLayers 3 with .net and the sample has this code for getting feature info from a WMS.

map.on('singleclick', function(evt) {
  document.getElementById('info').innerHTML = '';
  var viewResolution = /** @type {number} */ (view.getResolution());
  var url = wmsSource.getGetFeatureInfoUrl(
      evt.coordinate, viewResolution, 'EPSG:3857',
      {'INFO_FORMAT': 'text/html'});
  if (url) {
    document.getElementById('info').innerHTML =
        '<iframe seamless src="' + url + '"></iframe>';
  }
});

This line give me sintaxis error.

var viewResolution = /** @type {number} */ (view.getResolution());

But several sites have the same sample code.

I can see looks like a comment line.

So how i fix it.


texts: Filling Polygon with SVG File so it Will Appear Only One Time in QGIS 3.4.2

texts:
I have a layer that shows the world, without Antarctica and small islands.
This layer is made of five different polygons.

I want to fill it with an SVG file that is displayed non uniform background.

When I create a point layer the SVG background is displayed perfectly.

However, within the polygon I was unable to display it  properly.

However I change the Texture width parameters the SVG was not fitting perfectly.
What is needed to be done in order of make it fit?


texts: Downgrade to long term release QGIS with OSGeo4W?

texts:
QGIS this year split into two branches:

Latest release (now 2.10.2)

Long term release (now 2.8.2)

I am not sure how is it normally in larger companies but I think having (hard)core GIS employees who need to work with OSGeo4W version and GIS users who are fine with standalone version is not that rare. And if these two groups need to work on the same projects they need same versions of QGIS.

That makes the GIS employees need to have either two versions (OSGeo4W 2.10 and standalone 2.8) or possibility to be able to downgrade with OSGeo4W to long term release version.

Is downgrade with OSGeo4W to long term release possible somehow? At this stage I have possibility to only downgrade to 2.10.1 in the advanced installed list.


texts: Updating One Country's OSM Data through Osomosis

texts:
I would like to know how a single country's OSM data is kept up to date by osmosis. This is not a question regarding the method of doing the update, rather how osmosis does it.

I downloaded the OSM data for a single country and loaded it to a PostGIS database. Next, OSM minutely updates were setup. When I run the update, osmosis correctly updates only the data relevant to that specific country.

However, I can't figure out how it does that without any country or boundary information stored in the database.


texts: Calculating XY coordinates in Decimal Degrees from projected dataset in model or script?

texts:
Is it possible to calculate points XY coordinates in Decimal Degrees in model or script when a feature class is in projected coordinate system?

It's easy when a FC is in geographic coordinate system:

using Add XY Coordinates tool or
using Python expression eg. !shape.extent.XMax!

I've found that area and length properties of geometry field can be modified with geometry unit conversion keyword. For linear units of measue one can use @DECIMALDEGREES.
Unfortunatelly, !shape.extent.XMax@decimaldegrees! doesn't work as XMax is not a length.

In Calculate Geometry function (accessed from right-click) there's a possibility to choose Decimal Degrees output type even for projected feature class.
Can I do this using ArcPy?

Here's a code snippet based on iRfAn's solution:

import arcpy, os
projectedFC = r"C:\tmp\test.gdb\points01_Projected"
prjFile = os.path.join(arcpy.GetInstallInfo()["InstallDir"],
            r"Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj")
spatialRef = arcpy.SpatialReference(prjFile)

updCursor = arcpy.UpdateCursor(projectedFC,"", spatialRef)
for row in updCursor:
    pnt = row.Shape.getPart(0)
    row.X = pnt.X
    row.Y = pnt.Y
    updCursor.updateRow(row)

del updCursor, row


texts: gdal writeArray causes appcrash at end of python script

texts:
i'm using OSGeo4W64 (gdal 2.0.2-8) with Python 2.7.5 on Windows 8.1 64bit.
Editor is PyCharm 2016.1.2.

I want to save an image with gdal, but i always get an error when the python script is complete. The following example causes the error:

from osgeo import gdal
import numpy as np

path = 'D:\\test.tiff'
data = np.array([[1]])
sizeY = int(data.shape[0])
sizeX = int(data.shape[1])
bandCount = 1
driver = gdal.GetDriverByName("GTiff")

dataSet = driver.Create(path, sizeX, sizeY, bandCount, gdal.GDT_Byte)
band = dataSet.GetRasterBand(1)
band.WriteArray(data)  # causes error, if commented out, everything is fine

band = None
dataSet = None

print "the appcrash will come next"

I'm dereferencing the vars and it should be fine, but i'm getting "python.exe has stopped working".

The console shows:

the appcrash will come next

Process finished with exit code -1073741819 (0xC0000005)

When i run the code in QGIS (2.14.1) and i close QGIS, a similar error occurrs and a dump file was created.

Process Name:   qgis-bin.exe : C:\OSGeo4W64\bin\qgis-bin.exe
Process Architecture:   x64
Exception Code: 0xC0000005
Exception Information:  The thread tried to read from or write to a virtual
address for which it does not have the appropriate access.

My GDAL-Enviroment-Vars are:

GDAL_DRIVER_PATH=C:\OSGEO4~1\bin\gdalplugins
GDAL_DATA=C:\OSGEO4~1\share\gdal

I have no idea why this error happens (maybe something with the dereferencing). The error might be similar to here.


texts: QGIS cannot see rasters in postgres10 / postgis 2.4

texts: I am running a postgis-enabled postgres database. I recently did an upgrade from 9.6 to postgres10 using postgis 2.4 on both databases. all my vectordatasets still works fine, but I had a dtm as a raster data set in the data base that now shows up in without geometry. I have also tried to reimport the dtm, but to no avail. Is there anything special I need to do to activate the raster support? Raster calculations seems to work fine in postgres, so the problem must be somewhere in the interface between postgres and qgis. I am using qgis 2.14.11 and the nightly build of 2.99 - both with the same result.


texts: Multiple buffers to both sides of a line using ArcGIS

texts:
I am aiming to create multiple buffer zones to both sides of a line that have flat ends where the line ends (unlike 'rounded' edges of multiple ring buffers) and are not connected (each buffer zone to each side an individual feature) as shown for the blue line.

I am wondering if it is possible to do this without having to create each single buffer and then to clip and merge.


texts: Finding closest outside point (health center) to point (education center) inside polygon (danger zone) using PostGIS?

texts:
I am trying to find the closest health centers (points) to the education centers (points) located within the "danger zone" (a specific polygon) while excluding those health centers located in the "danger zone".

I'm pretty sure I need to join the three tables to do so...

so to get the education centers within the danger zone I have a virtual table:

   with edu_t as(
   select distinct on (e.gid) e.gid as edu_gid, e.name as edu_center, e.geom as edu_geom
   from education_centers e, danger_zones d
   where (st_contains(d.geom, e.geom)) and d.name = 'Danger Zone 1')

and to get the closest health center from these points I have:

   select distinct on (et.edu_gid) et.edu_gid, et.edu_center, h.gid, h.name,
   st_distance(et.edu_geom, h.geom) as dist
   from edu_t et, health_centers h
   where st_dwithin(et.edu_geom, h.geom, 50000)
   order by et.edu_gid, st_distance(et.edu_geom, h.geom) asc;

but that mainly gives me health centers in Danger Zone 1.

How can I exclude them?


texts: Assigning numbers manually in ArcGIS field calculator?

texts:
I want to assign 3 numbers manually in a field

at the field of Ddensity i want to assign 10,40,85 ... how can i do that?


texts:
Translation
When QGIS is installed, it displays the following message:

the ordinal 3109 could not be located in the dynamic library LIBEAY32.dll.

Could someone please help me.

Original
Cuando se instala QGIS, se despliega el siguiente mesaje:
"the ordinal 3109 could nor be located in the dynamic library LIBEAY32.dll."
por favor alguien me puede ayudar.


texts: Connected nodes and lines in a pipe network

texts:
I'm relatively new to GIS. I'm currently using QGIS.

I have 2 shapesfiles for a sanitary pipe network - a manhole file with the manhole nodes, and a pipe file with the pipe lines.

I need to find a way to identify the nodes which a line is "connected" to.

ie: Pipe ID n is connected to Manhole IDs x and y


texts: WMS layer option with Gemini VA WMS Server on Leaflet

texts:
I have got a WMS layer from an instance of Gemini VA WMS Server I want to add to my Leaflet map, version 0.7.7.
For that I need to break down my WMS link to fit the options like this example:

var wmsLayer = L.tileLayer.wms('http://demo.opengeo.org/geoserver/ows?', {
    layers: 'nasa:bluemarble'
}).addTo(map);

This is the link I am not sure how to use. In QGIS the first line is ID 0 and name "RootLayer". It contain 8 layers and all of them make up the layer I want to add.

The link I have provided is not live but it follows the same structure to the link I got that is working in QGIS.

EDIT: As suggested from comments I tried this:

var wmsLayer = L.tileLayer.wms('http://localhost/MapServer/Standard_SDO.wms', {
}).addTo(map);

That throws me an error: Uncaught TypeError: Cannot read property 'wms' of undefined.

It seem that it does not like the dot.
Do you know any work around?

This is what I have tried so far but it did not work:

var wmsLayer = L.tileLayer.wms('http://localhost/MapServer/ows?', {
        layers: 'Standard_SDO.wms'
    }).addTo(map);


texts: I am trying to calculate EVI from MODIS data in Arc. I know how to run one raster group of 3 bands but don't know how to tell it to move on to the next group of 3 bands in that directory. Here is the formula EVI = 2.5 * (NIR-Red)/(NIR+6*Red-7.5*Blue+1). So ideally I would like it to run a directory of many tiles with their bands in order from 1-3.


texts: Using ArcGIS Field Calculator to populate number field from text field?

texts:
I'm trying to assign speed values to street types using TIGER street shapefiles. So I created the new field for the speed, and I was going to use the field calculator to take the MTFCC classifications and fill the speed field from that, using a simple replace script.

I have this simple script that I found from another question that someone asked. I plugged my fields and values in and repeated the replace.

def Calc(value):
if "S110" in value:
    return value.replace('S1100','65')
if "S1200" in value:
    return value.replace('S1200','45')
if "1400" in value:
    return value.replace('S1400','30')
if "S1500" in value:
    return value.replace('S1500','15')
if "S1630" in value:
    return value.replace('S1630','25')
if "S1640" in value:
    return value.replace('S1640','25')

However, when I run it, it seems to run fine, but my fields are still Null. I really have no experience with Python, so I'm not even sure where I'm going wrong with the script.


texts:
I would like to hide/display the map when clicking somewhere. Does Folium have any parameters for that?
folium.Map(location=[48.86088, 2.33765], zoom_start=8)


texts: Acquire direction between points along route

texts:
How can I generate a directionality attribute between two points along a route/line?

Currently I am using Network Analyst to determine closest facilities, but I need a direction (upstream/downstream) along the route from incidents to facilities. Below is a visual example. The attribute I am trying to generate would be the direction ("Upstream") along the route you would have to travel to get from A to B. It should be noted that A and B are different features of different feature classes.

I am somewhat new to Network Analyst. Also, I am open to using methods other than Network Analyst.

EDIT:
Clarification


texts: Define raster origin so each instance of LiDAR interpolation has same grid

texts:
Is there a way to control the origin point of the resulting raster when interpolating from a LAS tile?

I have a large collection of regulate .las tiles with some overlap between tiles (1,000m tiling system with 50m overlap so 1,100m x 1,000m total tile size). I want to individually interpolate using IDW or Discrete Sibson as needed, but ensure that resulting rasters have the same underlying grid so that there is no need to re-interpolate when merging the rasters. As it stands, when i interpolate multiple files they all are on slightly different grids. My work involves tracing small artifacts like walls and ditches that get lost quickly if I process too much. Methods considered so far:

Create artificial LAS points at the exact corners of each original .las file and merge to create a file with the exact bounding box of the tiles
Do the same thing in Whitebox after converting each .las file to .shp

Both methods seem to work but require more automation then I have the skills for. If there is some way I am missing to just make each instance of LiDAR interpolation in Whitebox adhere to a grid that I set that would be much easier.


texts: Merge shapefiles into geodatabase in ArcPy

texts:
I am trying to merge a number of shapefiles with the following code:

import arcpy, os, datetime

#Set working directory
cd = r"C:\Users\..."
arcpy.env.workspace = cd
os.chdir(cd)
arcpy.env.overwriteOutput = True

directory = cd + "/shapefiles"
shapefiles = []
for root, dirs, files in os.walk(directory):
    files = [ fi for fi in files if fi.endswith(".shp") ]
    for filename in files:
        shapefile = (os.path.join(directory, filename))
        shapefiles.append (shapefile)

arcpy.Merge_management(shapefiles, cd + "all_districts.shp")

This fails due to Error 001156 which arises because some fields in some shapefiles exceed the maximum text length. According to this question, it is not possible to increase the maximum text length and the answers suggest that one should merge the shapefiles within a geodatabase.

I am new to geodatabases, so I do not know how to do this. How can I merge all shapefiles from my folder into a geodatabase?


texts: Qgis server - tiles are not rendered when version is set to 1.3.0

texts:
I am playing with qgis server and just plain openlayers 3 as a client. With the following code i get the map ro render properly:

var projection = new ol.proj.Projection({
      code: 'EPSG:32633',
      units: 'm',
      axisOrientation: 'neu'
  });

var map = new ol.Map({
        target : "map-panel",
        layers : [
          new ol.layer.Tile({
              source: new ol.source.TileWMS({
              url: 'http://192.168.200.142:8081/cgi-bin/project/qgis_mapserv.fcgi',
              params: {'FORMAT': "image/png",
                       'VERSION': '1.1.1',
                       tiled: true,
                       LAYERS: 'roads',
                       STYLES: '',
              }
            })
          })
        ],
        view: new ol.View({
          projection: projection
        })
    });

If i change the VERSION parameter of the WMS request to 1.3.0, qgis server returns blank tiles (requests are processed). Why is that? I am running 2.14.0 64bit on Linux Mint

EDIT: Here is the actual request URL (when version is set to 1.3.0):

http://192.168.200.142:8081/cgi-bin/project/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&tiled=true&LAYERS=roads&WIDTH=256&HEIGHT=256&CRS=EPSG%3A32633&STYLES=&BBOX=5472872.718699152%2C625471.1678513326%2C5629240.510661985%2C781838.9598141655' -H 'Accept: image/png,image/;q=0.8,/*;q=0.5' -H 'Accept-Encoding: gzip, deflate' -H 'Accept-Language: en-US,en;q=0.5' -H 'Connection: keep-alive' -H 'Cookie: session_id=ac18ac6c07b8da6034fae3c3ae6a5a3ed7e61f87' -H 'Host: 192.168.200.142:8081' -H 'Referer: http://192.168.200.142/apps/app' -H 'User-Agent: Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:45.0) Gecko/20100101 Firefox/45.0'


texts: How to style virtual vertices when using the ModifyFeature control with OpenLayers?

texts:
I am building an Android application that uses OpenLayers embedded in a WebView.  Being mindful of the imprecision of touch interfaces I am trying to style all controls to be large and easy to operate.

I have followed the documentation on feature styling but have been unable to control the size of the virtual vertices when using the createVertices mode of the ModifyFeature control.

By way of example consider the following style with a pointRadius of 5:

var customStyle = new OpenLayers.Style({
                    fillColor: "black",
                    strokeColor: "red",
                    pointRadius: 5,
                    strokeWidth: 30
                });

var customStyles = new OpenLayers.StyleMap({
                    "default": customStyle,
                    "select": customStyle,
                    "temporary": customStyle
                });

var polyLayer = new OpenLayers.Layer.Vector("Polygon_Layer", {styleMap: customStyles });

Meanwhile the following has a pointRadius of 50:

var customStyle = new OpenLayers.Style({
                    fillColor: "black",
                    strokeColor: "red",
                    pointRadius: 50,
                    strokeWidth: 30
                });

As you can see this has no impact on the size of the virtual vertices midway along each edge.

Does anyone know a technique to control their size?


texts: Using project on point shape file to change coordinate system using ModelBuilder

texts:
I have created a model using ModelBuilder in ArcMap 10.8 to do the following:

Project point shape file to change coordinate system from A to B
Add a field (X)
Calculate field: !SHAPE.CENTROID.X! using python_9.3
Add a field (Y)
Calculate field: !SHAPE.CENTROID.Y! using python_9.3
Table to excel.

Whenever I run the model it stops at "Project" and gives me the following error ERROR 999999: Error executing function. Create output feature class failed Failed to execute (Project)."
Why does it not work? When I do this manually it works.


texts:
I periodically need to refresh my QGIS map based on changed underlying data sources (currently locally-saved KML, GPKG, but trying to be general) using PyQGIS. Following Make a Qgis layer update from a changed data source, for a specific layer I can use:

layer.dataProvider().forceReload()
layer.triggerRepaint()
iface.layerTreeView().refreshLayerSymbology(layer.id())

However, I need to find all current layers depending on a specific updated source, and also keep track of the layers found vs the layers (=database tables, KML folders) in the source. The latter is to alert the user if there are map layers to deprecate that are no longer in the datasource, or new ones in the datasource to potentially add and style. For this I need to parse layer.source() to extract filename and (source as opposed to map) layer name.

What's a robust, non-provider specific, and reasonably future-proof way of doing this? At the moment, all the relevant sources are loaded via ogr, and by inspection the relevant sourcestring has the URI-like format full_filename|layername=table_or_folder_name_in_source and sometimes also  ...|geometrytype=Point_or_LineString_or_Whatever. Therefore I use a simple if filename in layer.source() to find layers that need updating, and layer.source().split('|layername=')[1].split('|')[0] to pick out the layer name in the source file. But this could be confused by layers with different providers (e.g. a WMS layer whose URI happens to have a confounding substring) and for all I know ogr could change its datasource string format at any time.

More specifically, is there either a method/function to parse a layer.source() string into meaningful components, or -- failing that -- a documented canonical format for the string? I've noticed some providers use & as a delimiter but others |.


texts: Make a line chart using the ui.Chart.array.values function (Earth Engine)

texts:
Is there a way to make a line chart with the ui.Chart.array.values function? When I try to set the chart type using .setChartType('LineChart'), it shows only points and has a weird scale on the h-axis.

Additionally, I want to include another band for a different year, say VIIRS2019 to show two plots on the same chart and do a visual comparison, can I just add the 'avg_rad' from that other image into my VIIRS2014?

// CITY AND TOWNS SHAPEFILES
var point =
    /* color: #d63000 */
    /* shown: false */
    ee.Geometry.Point([-118.243683, 34.052235]);
var cities = ee.FeatureCollection("users/j9herrer/USACities2");

var mycity = cities.filterBounds(point);
//var la = ee.Geometry.Point([-118.243683, 34.052235]);
mycity = mycity.filterBounds(point);
var city_outline = ee.Image().byte().paint({
  featureCollection: mycity,
  width: 2
});

Map.addLayer(mycity);

var latLon = ee.Image.pixelLonLat().select('longitude');

// TRANSECT
var geom_city = mycity.geometry();

var boundingbox=geom_city.bounds();
var coords = boundingbox.coordinates().get(0);
//for whatever reason .coordinates returns a list with one element which is the list we want, hence .get(0)
coords = ee.List(coords);
var NW = coords.get(1);
var SE = coords.get(3);
var mydiagonal = ee.Geometry.LineString([NW,SE]);
Map.addLayer(mydiagonal);
print (mydiagonal, "Transect");

// VIIRS DATA
var VIIRS = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG");
VIIRS = VIIRS.select('avg_rad');

var VIIRS2014 = VIIRS.filterDate('2014-1-1','2015-1-1').mean().addBands(latLon);
var VIIRS2019 = VIIRS.filterDate('2019-1-1','2020-1-1').mean().addBands(latLon);

var VIIR2014List = VIIRS2014.reduceRegion({
  reducer: ee.Reducer.toList(),
  geometry: mydiagonal,
  scale: 1000,
  maxPixels: 1e9
});

print (VIIR2014List, "VIIRS 2014, List");

var lon = ee.List(VIIR2014List.get('longitude'));
var avgRad = ee.List(VIIR2014List.get('avg_rad'));
var chart = ui.Chart.array.values(avgRad, 0, lon ).setChartType('LineChart');

print(chart);


texts:
I'm having a problem with the following script.  When listing parameters that are optional I've used the "" before, so I can't figure out why this isn't working.

import arcpy
import os

arcpy.env.overwriteOutput = True

path = "C:/Data_new/Temp/Default.gdb"
name = "point_shapefile.shp"
type = "POINT"

m = "DISABLED"
z = "DISABLED"
prjFile = "C:/Program Files/ArcGIS/Desktop10.0/Coordinate Systems/Projected Coordinate Systems" + \
      "/State Plane/NAD 1983 (US Feet)/NAD 1983 StatePlane North Carolina FIPS 3200 (US Feet).prj"
reference = arcpy.SpatialReference(prjFile)

arcpy.CreateFeatureclass_management(path, name, type, "", m, z, reference)


texts: Creating raster image using Shapefile in ENVI?

texts:
I have shapefile containing polygon boundaries and want to change them as a raster image in ENVI 5.3. Although, this could easily be done in ArcMap or other GIS software, however, I will have to perform all operations in ENVI only.

I have searched ENVI toolbox to see if there is any availability but didn't succeed.


texts: GPS Elephant Tracking: Utilizing A-GPS, WAAS, DGPS and Understanding HDOP, NSAT, and ZDA

texts:
*I'm using a GPS data logger (qstarz bt-q1000xt) to track the movement of zoo elephants, and I have a variety of technical questions.  I've had limited success receiving answers from the vendor.  You can see the GPS anklet on an elephant (links to flickr) here.

A-GPS:
In mobile phones, A-GPS refers to the ability to use cell towers to triangulate position when SV's are not visible.  In GPS units, A-GPS means you can download almanac data via the internet for faster TTFF, rather than downloading the almanac from a satellite.  Is this a correct description?  If so, does this mean that A-GPS provides only improved TIFF, and has no effect on coordinate accuracy?  The reason I ask is that I want to send these units out to other zoos to collect data, but don't want them to have to screw around with the software.

WAAS/DGPS:
There seems to be varying opinions about whether or not WAAS is a type of DGPS.  For example, my data loggers state they have "DGPS(WAAS+EGNOS+MSAS) support."  Yet my understanding is that WAAS is an SBAS that my device can access and utilize in real time, and DGPS is an GBAS which is strictly a post-processing analysis done with software.  But perhaps this is wrong, and you can purchase fancy receivers that allow your unit to execute differential correction in real time?  Even if this is true, for budgetary reasons I am more interested in post-processing with DGPS.  What level of expertise is required to conduct DGPS analysis?  Where would be a good place to start looking for more information on this?  The reason I ask is that I hope to evaluate elephant herd dynamics, for which I need meter or sub-meter accuracy, which my data logger alone cannot seem to provide.

HDOP/NSAT:
Zoo researchers using GPS regularly remove fixes with an HDOP less than 2, and will then also remove fixes with an NSAT less than 6, "just to be safe."  However, 4 satellites with perfect geometry should be able to make a highly accurate fix.  Is there any reason to remove fixes with NSATS of 4 or 5?  I also understand that an HDOP value of 1 is "perfect;" meaning the satellite geometry is multiplying the effect of other errors (atmospheric effects, satellite clock error, etc) by a factor of 1 (i.e., not at all).  But my data loggers occasionally give an HDOP value of less than one.  How can this value be interpreted?

ZDA:
Are the times provided by this NMEA sentence more accurate somehow than the times given by RMC?  Why would you want to record ZDA?

An associate told me that the irony of GPS is that everyone knows enough to be dangerous, but nobody REALLY knows how it works.  Prove him wrong?

Thanks!
Matty


texts:
Possible Duplicate:
  How to create a point along a line given distance
  How to create equidistant points in QGIS?

I have some road survey data and each line has an individual value and defect data (which is given as a figure in meters of how far down that line it is) and I was wondering if it would be possible to generate points at the chainages given for these defects in relation to the lines they are on using the direction the line has been drawn in.


texts: Georeferencing an image with GDAL, based on QGIS Georeferencer Plugin

texts:
I would like to georeference an image with GDAL. The image already has a spatial reference set (EPSG:31276). For the georeferencing, I have a list of 4 ground control points (also in EPSG:31276).

Using the QGIS Georeferencer Plugin, I can correctly georeference the image. I then select "Generate GDAL Script", which results in the following GDAL commands:

gdal_translate -of GTiff -gcp 6.55253e+06 -5.07402e+06 6.55273e+06 5.07395e+06 -gcp 6.55246e+06 -5.0711e+06 6.55263e+06 5.07106e+06 -gcp 6.55477e+06 -5.07382e+06 6.55511e+06 5.07368e+06 -gcp 6.55483e+06 -5.07197e+06 6.55525e+06 5.07197e+06 "C:/gis_sample_data/sample_image/image1.tif" "C:/Users/XYZ/AppData/Local/Temp/image1.tif"

gdalwarp -r near -order 1 -co COMPRESS=NONE  "C:/Users/XYZ/AppData/Local/Temp/image1.tif" "C:/gis_sample_data/sample_image/image1_georef.tif"

Those commands run without an error, but the resulting image is not correctly projected. Here is a small part from gdalinfo "C:/gis_sample_data/sample_image/image1.tif" (the original image):

Corner Coordinates:
Upper Left  ( 6551897.000, 5074277.000) ( 18d40' 4.21"E, 45d48'39.73"N)

whereas gdalinfo "C:/gis_sample_data/sample_image/image1_georef.tif" (result from the GDAL call) shows

Corner Coordinates:
Upper Left  ( -534350.794,  156938.268) ( 34d51'24.74"W,  0d52'47.63"N)

How would I have to modify the GDAL commands to correctly georeference the image?

The sample image is part of the Baranja Hill data set, which can be downloaded from here.
The ground control points are those:

mapX,mapY,pixelX,pixelY,enable
6552733.64613180514425039,5073946.04477077350020409,6552533.84482758585363626,5074021.48275862075388432,1
6552633.32449856773018837,5071059.01110315229743719,6552461.44827586226165295,5071104.32758620660752058,1
6555111.64040114637464285,5073678.52041547279804945,6554765.36206896603107452,5073821.32758620660752058,1
6555249.11819484271109104,5071965.62141833826899529,6554833.5,5071973.08620689623057842,1

When copied into a text editor and saved with the .points extension, they can be imported in the QGIS Georeference Plugin.

EDIT: Changed negative to positive values in gdal_translate. Results of gdalinfo:

Corner Coordinates:
Upper Left  ( -591887.923,  223739.838) ( 35d 8'28.93"W,  1d14'58.74"N)

By the way, I also tried adding -t_srs EPSG:31276, -s_srs EPSG:31276 and -a_srs EPSG:31276 and various combinations thereof to the gdal commands, but it didn't fix the problem.


texts:
First: Grass did not run right out of the box when installed as instructed through the QGIS 3.0.x package. I had to install Grass separately and create a database path.

Now I am getting the following error: tput: No value for $TERM and no -T specified

Thoughts?

See log message below:

    Welcome to GRASS GIS 7.2.2

    GRASS GIS homepage: http://grass.osgeo.org

    This version running through: Bash Shell (/bin/bash)

    Help is available with the command: g.manual -i

    See the licence terms with: g.version -c

    See citation options with: g.version -x

    Start the GUI with: g.gui text

    When ready to quit enter: exit

    tput: No value for $TERM and no -T specified

    Cleaning up temporary files...

    Done.

  Goodbye from GRASS GIS

  Default locale not found, using UTF-8

  Execution completed in 2.31 seconds Results: {'output':
  }

  Loading resulting layers The following layers were not correctly
  generated.
  /var/folders/tm/6kzq8c3d28z9y85wf0bw579h0000gn/T/processing_7475f0c4cfc745e981e77acc2a1ae6fa/464bb56c81354bf89fefbf6531a83898/output.tif
  You can check the 'Log Messages Panel' in QGIS main window to find
  more information about the execution of the algorithm.


texts: Performing spatial join in Spatialite (QspatiaLite for QGIS)

texts:
I'm trying to spatially join and summarize the attributes of a set of buildings to a set of containing parcels. There are (potentially) many buildings to each one parcel.

Here is the query as I'm attempting to run it:

SELECT "MF_with_Zones_atts".'Geometry' AS Geometry,
"MF_with_Zones_atts".'PKUID' AS ID,
"Apt_Bldg_Footprints".'STATE_ID' AS STATE_ID,
"MF_with_Zones_atts".'JURIS_CITY' AS JURISDICTION,
"MF_with_Zones_atts".'MAIL_CITY' AS CITY,
"MF_with_Zones_atts".'STATE' AS STATE,
"MF_with_Zones_atts".'ZIPCODE' AS ZIP,
"MF_with_Zones_atts".'COUNTY' AS COUNTY,
"MF_with_Zones_atts".'UNITS' AS UNITS,
"MF_with_Zones_atts".'UNIT_TYPE' AS UNIT_TYPE,
"MF_with_Zones_atts".'ZONE' AS ZONE,
"MF_with_Zones_atts".'CMP_' AS CMP,
"MF_with_Zones_atts".'OVRLY' AS OVERLAY,
"Apt_Bldg_Footprints".'BLDG_USE' as BLDG_USE,
"MF_with_Zones_atts".'ZoneMaxDU' as ZoneMaxDU,
"MF_with_Zones_atts".'ZoneMaxFAR' as ZoneMaxFAR,

count("Apt_Bldg_Footprints".'ROWID' AS BLDG_COUNT,

sum("Apt_Bldg_Footprints".'BLDG_SQFT') as BLDG_SQFT,
sum("Apt_Bldg_Footprints".'UNITS_RES') as BLDG_UNITS_RES,
sum("Apt_Bldg_Footprints".'Bldg_Ftprt') as BLDG_FOOTPRINT,
max("Apt_Bldg_Footprints".'NUM_STORY') as NUM_STORY,
max("Apt_Bldg_Footprints".'AVG_HEIGHT') as AVG_HEIGHT,
max("Apt_Bldg_Footprints".'MAX_HEIGHT') as MAX_HEIGHT,

FROM "MF_with_Zones_atts", "Apt_Bldg_Footprints"
WHERE Contains("MF_with_Zones_atts", "Apt_Bldg_Footprints")
GROUP BY "MF_with_Zones_atts".'ADDRESS'
ORDER BY "MF_with_Zones_atts".'ADDRESS'

The error messages I get are:

First:

  "The SQL query seems to be invalid. near "AS": syntax error

Then:

  Unable to read temp table: MF_w_Footprints.tmp Operation cancelled

So, the question is, what am I doing wrong?

Both layers have been imported into a new sqlite database as geometry tables, and I'm using the QspatiaLite plugin for QGis to try to make this happen.


texts: How to convert an image of a floor plan into data that I can work with?

texts:
This is more of a research question. I'm looking for aways of converting an image of a floor plant into data that I can work somehow. The end goal is to something similar to this: http://indoorosm.uni-hd.de/3d/.

I want something that allows me to store rooms and corridors informations, and then use it to trace routes between them.

It would be really helpful is someone can point me in the right direction, because the more I search the more different technologies I find, and I don't know which on would serve best my goals.


texts: Getting feature from coordinates in geotools

texts:
I'm working on web services in Java and using GeoTools java api for wfs operations.

I have a wfs store and a lot of features in it. I do some operations via geotools wfs api for example buffer, intersect etc. My problem is not building feature from points. I'm using geotools for building my geo web services and want to create a web service method that must use two parameters(lat and lon) and return content of a feature.

This returned feature must be perfectly matched with given latitude and longitude so this coordinates should point out this feature

How can I do that ?

Such as wms,wfs ogc services are rest apis and return xml responses. I don't want to use this services  with http request or rest api clients.Also I don't want to parse directly response document at the same time doing all operations in client side is not useful for my app. I want to build my web service composition for some defined operations in wfs,wms and etc.


texts: I am trying to merge 166 tiles into a single world raster (500 meters resolution) through ArcGis "Mosaic to new raster" tool with a simple Python script. The process seems incredibly slow. Before starting to study alternative approaches (from similar questions GDAL_merge seems the most correct) is there any ArcGis setting that can speed up the process?


texts: Unable to create .esriAddin; missing ESRI ArcGIS Add-in SDK components(s)?

texts:
What am I missing?  I upgraded to ArcGIS 10.3 succesfully recompiled all other 20 or so addins after encountering a similar issue and renaming "ESRI.ArcGIS.AddIns.11.targets" to "ESRI.ArcGIS.AddIns.targets" possibly due to a known defect (#01599069)

Now I would like to create a brand new addin but I get the same warning. If if rename "ESRI.ArcGIS.AddIns.targets" back to "ESRI.ArcGIS.AddIns.11.targets" the warning dissapears and addin seems to compile but none of my old addins will recompile raising the same old error.

If i copy and paste the "ESRI.ArcGIS.AddIns.targets" and also keep the "ESRI.ArcGIS.AddIns.11.targets" both the old and new addins compile without errors or warning but the new addins ref ID are not recognized by the older upgraded addins.

I am pretty much completely confused as to what needs to be done here to make this work.