texts: Loading Postgres layer with Session parameter (SET) using Pyqgis

texts:
With postgreSQL 9.5 and +, there is a policy feature that handle row-level access which is great.

My special case is to use this feature to control access to a multi-user postgres DB.

Here is what i'm trying to do:

I have a postgres table (let's say test_data) with a policy associated. The policy is defined as follows:

CREATE POLICY test_data_pol ON test_data FOR ALL USING ('conc' = get_pol())

conc is a column in test_data which contains integer representing id of user that have access to this row of data.
get_pol() is a postgres function that retrieve the value of a session parameter (my.foo) set by another function set_pol() as follow:

set_pol()

CREATE OR REPLACE FUNCTION public.set_pol(foo_val integer)
  RETURNS void AS
$BODY$
BEGIN
  EXECUTE format('SET SESSION my.foo TO %I', foo_val) ;
END;
$BODY$

get_pol()

CREATE OR REPLACE FUNCTION public.get_pol()
  RETURNS character varying AS
$BODY$
DECLARE
  foo_val character varying;
BEGIN
  SHOW my.foo INTO foo_val;
  RETURN foo_val;
END;
$BODY$

With that system I can handle the access to the row by executing a SELECT set_pol(value) when I connect to the DB. After setting this Session parameter I can access only the relevant rows.

This works well with psycopg2. for example:

import psycopg2
conn = psycopg2.connect(host=DATABASE['host'],
                        database=DATABASE['dbname'],
                        user=DATABASE['user'],
                        password=DATABASE['password'],
                        port=DATABASE['port'])
cursor = conn.cursor()
cursor.execute("SELECT set_pol(%s)" % conc_value)
cursor.execute("""SELECT * from test_data""")
data = cursor.fetchall()

With that python snippet data only contains the row that are accessible with the relevant conc_value.

Now I try to add VectorLayer in QGIS using pyqgis with this code snippet:

uri = QgsDataSourceURI()
uri.setConnection(host,
                   port,
                   dbname,
                   user,
                   pasword)
uri.setParam('my.foo', 'conc_value')
uri.setDataSource('public', 'test_data', None, None, 'id')
lyr = QgsVectorLayer(uri.uri(), 'test_data', 'postgres')

But the setParam statement does not create the variable as expected in the connection Session and I get this error message:

ERROR:  unrecognized configuration parameter "my.foo"
CONTEXT:  SQL statement "SHOW my.foo"
PL/pgSQL function get_pol() line 5 at SQL statement
.
SQL : SELECT * FROM "public"."test_data" LIMIT 1

which make me think that the setParam function does not do what I need... I need to execute an SET SESSION my.foo TO 'conc_value' but setParam doesn't seems to  create the variable the proper way as SHOW my.foo in get_pol() raise an error that it does not know the my.foo parameter...

So my question is:

How can I create the connection with pyqgis to handle the creation of my.foo parameter at the session level?

OR

How can I modify my get_pol() function to be able to access to the my.foo parameter create by the QgsDataSourceURI.setParam() pyqgis function?


texts: Image Dodging (color matching algorithm) in ERDAS

texts:
I have two LISS 3 images acquired in Feb and March a year apart. I have performed atmospheric correction (haze reduction) for both the images individually before applying the color matching algorithm while mosaicing. I applied the image dodging color matching algorithm in the mosaicpro tool. The output was fine, except that the edges of the two images seemed to be lightened, i.e., the mosaic was not seamless in the edges. It could be a radiometric error for the edges of the two images. I could not find any specific parameters in the tool to overcome this error.
I have tried this procedure (color matching two images before mosaicing) in ArcGIS and ENVI but the results were not satisfactory.
Can I overcome this error using ERDAS?


texts: Unable to load Vector Layer in standalone application developed using QGIS libraries?

texts:
I suspect that the problem is due to the Prefix path. My ubuntu system does not have it set. What to do?

>>> print(QgsApplication.showSettings())
Application state:
QGIS_PREFIX_PATH env var:
Prefix:
Plugin Path:
Package Data Path:
Active Theme Name:
Active Theme Path:  /themes//icons/
Default Theme Path: :/images/themes/default/
SVG Search Paths:
User DB Path:   /resources/qgis.db
Auth DB Path:   qgis-auth.db


texts: Automating recording fields in a shapefile via pyshp

texts:
I am following the tutorial here for how to convert a csv to shp using pyshp. I am almost there, with one issue. In the tutorial, the person loops over the rows in the csv, and for each field, assigns the value of that field to a variable (shortened example):
import shapefile # import pyshp library
trees_shp = shapefile.Writer(shapefile.POINT)
for row in reader:
   tree_id = row[0]
   address = row[1]
   town = row[2]
   ....
   condition = row[11]
   latitude = row[12]
   longitude = row[13]

They then record the new point with attributes like this:
trees_shp.point(float(longitude), float(latitude))
trees_shp.record(tree_id, address, town, ..., condition)

My question is, is there another (better) way to record field values that are stored in a list (or comma separated string)? i.e. I am looking to record the attributes that are stored in 2 variables (for each row/point):
row_str = 'fieldValue1, fieldValue2, fieldValue3'
row_list = ['fieldValue1', 'fieldValue2', 'fieldValue3']

I would like to just do something like:
outShpWriter = shapefile.Writer(shapefile.POINT)
outShpWriter.record(row_str)
# or
outShpWriter.record(row_list)

but each gives me the following error:
  File "/usr/local/lib/python2.7/dist-packages/shapefile.py", line 1071, in record
    record = [recordList[i] for i in range(fieldCount)]
IndexError: tuple index out of range

Is there a workaround besides having to explicitly define variables i.e. (f1 = row_list[0]) and record those variables to the shp?
I ask because the number of fields/fields themselves may change and I would like to automate as much as possible to make it flexible.


texts: QGIS 3.6 splitwithlines algorithm does not always split

texts:
I'm trying to split one line layer with another line layer, and additionally get the vertices where they intersect.

However, if a line is self-intersecting, it does not split the underlying layer. In addition, some of the lines that are touching the underlying line, but not crossing them are not creating a split. Check the gist with the data to try it yourself.

I guess some floating point precision problem here, so I'm open for any advice on how to solve this issue.

As a workaround, I've tried extending the lines, but then I get two intersection vertices where the two lines meet on a baseline.

Check the sources here: https://gist.github.com/suricactus/9e43005bc1ca5f35d375bbde787a8a8e


texts: How to insert layer excluding MXD files using arcpy?

texts:
I'm working on 24 MXD files. Several maps include "roads" layer in them, the other maps don't have this layer. I'm trying to add "roads" layer to all the maps that do not have this layer by using this code:

import arcpy,os,sys,fnmatch
from arcpy import env

env.workspace = r"G:\desktop\Project"
insertLayer = arcpy.mapping.Layer(r"G:\desktop\Project\layers\roads.lyr")
counter = 0

for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"G:\desktop\Project\\" + mxdname)
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    refLayer = arcpy.mapping.ListLayers(mxd, "residence", df)[0]
for lyr in arcpy.mapping.ListLayers(mxd, "" ,df):
    if 'roads' not in lyr:
        arcpy.mapping.InsertLayer(df, refLayer, insertLayer, "BEFORE") # BEFORE\ AFTER
mxd.save()
counter = counter + 1

del mxd
print counter

and this is the error i get:

 >>>
landUse.mxd
Traceback (most recent call last):
  File "G:/desktop/y/yaron/shonot/software/---gis---/tools/YARON_SCRIPTS/aaaaaaaaaaaaaaaaa.py", line 18, in <module>
if 'roads' not in lyr:
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\mixins.py", line 361, in __iter__
self.layers
AttributeError: LayerObject: Get attribute layers does not exist
>>>

UPDATE:

I added this code line:

  if 'roads' not in lyr.name:

and now 'roads' layer has been added several times (each map) to all maps no matter if they contain 'roads' layer or not.


texts: How do I perform a 3 parameter geotransformation and projection in Sql Server 2012?

texts:
I have a table with latitude longitude (NAD27) columns in it.  I compute two other columns, X and Y, representing Web Mercator (WGS84) location.

Currently I'm using a Arcmap to do this, by applying the recommended geotransformation for the study area - the 3 parameter (geocentric) geotransformation - to go from NAD27 to WGS84.

I would like to do this entirely within Sql Server 2012. From what I can tell, Sql Server does not support datum transformations out of the box. Does anyone know of a Sql library that supports this geotransformation?  I would like to simply use the same coefficients in Sql that I'm currently using in Arcmap.

I also need to project from WGS84 lat/long into web mercator.  I see this formula implemented in javascript, but if someone has a Sql stored procedure that does this, it would be great.


texts: Fill in areas between lines with polygons?

texts:
I have to make a geographical map.

I have to fill in the gaps between the lines with polygons. Is there a possibility to just trace those lines automatically?

I have to do everything manually now, which takes a lot of time.
Right now I am using the "Add object polygon" from the digitalizing toolbar.


texts: Get EditText value with event handler in MapBasic

texts:
I have a dialog in MapBasic that includes an EditText control, a PopupMenu control and a button. When the user presses the button a sub is called. In this sub I would like to get the values of the EditText and PopupMenu, do some maths with them and then output the results into a second EditText control (which will be disabled so as to be used for output only).

My question is how do I get the values from the EditText and Popupmenu controls?


texts: Printing Raster extents but not showing on console

texts:
print("Data {0}, xmin{1}, ymin{2}, xmax{3}, ymax{4}".format(ext.Xmin, ext.Xmax, ext.Ymin, ext.Ymax))

The line of code above is not printing out the extent for each raster in my dataset in the code below. It runs without an error, but nothing is displayed.

myRasterList = []
for RAWR in myRasterList:
    ext = arcpy.Describe(RAWR).Extent
    print(ext.Xmin)
    print ("Data extents {4}, Xmin{0}, Xmax{1}, Ymin{2}, Ymax{3}".format(ext.XMin, ext.XMax, ext.YMin, ext.YMax, RAWR))

import arcpy, sys, os
env = arcpy.env
from arcpy import management as DM
from arcpy.sa import *
from arcpy import analysis as AN
from arcpy import conversion as CV

env.workspace = r"F:/Fall 2019/GEOG 625/Exam Data/ExamData/MichaelBreidert/Pen/HMB.gdb"

scriptdir = os.path.dirname(sys.argv[0])
#ws = arcpy.env.workspace = r"F:\Fall 2019\GEOG 625\Exam Data\ExamData\MichaelBreidert"
# creating variable to hold GDB
hmb = scriptdir + "/HMB.gdb"

myRasterList = []
for RAWR in myRasterList:
    ext = arcpy.Describe(RAWR).Extent
    print("Data {0}, xmin{1}, ymin{2}, xmax{3}, ymax{4}".format(ext.Xmin, ext.Xmax, ext.Ymin, ext.Ymax))

How do I print the extents of all data in my HMB.gdb?


texts:
I have about 1000 vehicle equiped by a GPS sensor and 3G connection.
Each vehicle feeds a csv file  -using 3G connection- with its current location every 5 sec.
Each vehicle have its own CSV file.

At the end of the day, I end up with about 1000 CSV file, each one contain the vehicle track of the day.

My aim to export all the csv files to a map in a way to see the most frequented routes by time. For example at 5PM I want to have map showing the different roads colored from red to green to show the trafic flow rate at this time.

Of course the coordinates are not always so precise, so I will take into considerationerror margin. I will assume that any point outlying by 10m of a road is on that road.

Also, I have the possibility to create a database with all the 1000 vehicles data on it.

Now I think you have the big picture of the situation, if anyone could help and suggest a tool to do this.


texts:
I have a polygon layer. I extracted the coordinates of vertices, and now I need to add these coordinates in each polygon. The final shapefile should have an attribute table as follows:

where X1, X2, etc.. are the vertices coordinates for each polygon.

Any ideas on how to do this using ArcGIS for Desktop?


texts: Coordinate system disappears when I use publish map and iframe tag in GeoExplorer?

texts:
I managed to add xy coordinate in GeoExplorer, I use OpenGeo suite to use GeoExplorer. I added the OpenLayers mouse position code to the composer in order to show the coordinates. Please consider my previous question and answer:

how to show coordinate on a geoexplorer map by moving the mouse?

The problem is, when I publish my map and use iframe tag the coordinate system will disappear. I even changed the "viewer" to "composer" in iframe tag, but still no result.


texts: Export Polygons from Service Area tool using ModelBuilder?

texts:
I have created a service area in ModelBuilder.

I have assigned a polygon to the "facilities".

How do I export the service area polygon to a shapefile in ModelBuilder?


texts: Listing all Feature Classes being used in ArcGIS Pro project using ArcPy

texts:
I am trying to work out how to list all feature classes in ArcGIS Pro and ArcGIS Pro files. I have used this script in the past which works for .mxd but recently I have started using ArcGIS Pro which has another extension file. Can anyone clarify if they have attempted this in ArcGIS Pro as I have over 1000 feature classes that I would like a list of that are being used.
import arcpy
import os

mxdDir = r'C:\path\to\directory\of\MXDs'

mxdDict = {}

for mxd in os.listdir(mxdDir):
    mxdPath = os.path.join(mxdDir, mxd)
    mxdDoc = arcpy.mapping.MapDocument(mxdPath)

    for layer in arcpy.mapping.ListLayers(mxdDoc, "", ""):
        if layer.supports("DATASOURCE"):
            ds = layer.dataSource
            if ds not in mxdDict:
                mxdDict[ds] = [mxd]
            elif mxd in mxdDict[ds]:
                pass
            else:
                mxdDict[ds].append(mxd)


texts: Finding output of geoprocessing operation in ArcGIS for Desktop?

texts:
I am merging some features using the Geoprocessing Merge tool in ArcGIS 10.1 for Desktop. The tool requires that the output is saved in a folder. Although this was done after the process ended successfully i could not locate the final output in the folder i saved it in. I have tried searching all the other folders but still cannot find it.

How do i locate the final output?


texts: Error Message: The coordinates or measures are out of bounds

texts:
When extruding between 2 TINs using a polygon feature class, above error is received.

The coordinate system for each layer is consistent, the shapefile has no Z or M values, and I have recalculated XY extent.

But the error is still being returned. Is there anything I could be missing?


texts:
When you load an image service to ArcServer, you have the option to place an Item Description with the service.  This acts as the metadata, however it displays a duplicate: "Service Description"  then right below it "Description"

The information is an exact replication.

When you edit the Item Description of the mosaic dataset and restart the service, both Description blocks update.

How do I edit/remove one of the Description blocks to avoid redundancy?


texts: Designing interactive map of botanical garden in QGIS?

texts:
Using QGIS is absolutely new to me but I hope, we can make it as short and easy as possible to put my idea into practice. I hope that QGIS is the appropriate program for this purpose.

I have COREL-drawn (Coral Draw X6) map of a botanical garden, saved as SVG, which I need for my research. My intention is, to build up an interactive map/database, incorporating all my important species and their localities in the garden. And I was said, that the whole request is too long, so I'll start with the first step:

How can I import my map (as .svg or another format), that it keeps looking like it does in Corel? I found some tutorial which says to export the map as .DXF from Corel and then import this into QGIS, but the result looks pritty weird.


texts: absolute paths don't work for my leaflet custom icons

texts:
I've made an example leaflet project using a custom icon loading an image from the internet:

https://jsfiddle.net/dvdjsph/v30vL109/

As soon as I try to load something from my hard drive, however, the image won't load.  This is how I'm trying to load it.

var offerIcon = new OfferIcon(
    { iconUrl: "C:\image.jpg"});

L.marker([0, 0], {icon: offerIcon}).addTo(markers).bindPopup(name);

I know that JS fiddle won't load this -- it's on my hard drive.  I just put up the example to show that, minus that path, everything else is working fine.

Is there something I'm missing about paths with Leaflet?

EDIT:

At Matej's suggestion, I tried putting an icon in the same folder as my template and referencing it like so:

var offerIcon = new OfferIcon(
    { iconUrl: "image.jpg"});

Still not working.


texts: Which tool or plugin for Theissen / Voronoi polygons in QGIS?

texts:
As I understand you can generate Voronoi polygons in QGIS (2.18.2) using the following approaches;

Vector menu / Geometry tools /
OR
MMGIS plugin in

I also read here that you could use SAGA as well. I'm trying to create a shapefile of the UK postcode units using the Voronoi method as discussed here and wanted to understand what the difference was between the different algorithms given that I'll be running it across ~1.7 million points.


texts: Permanent override of webgis HTML file with Chrome: is it really possible?

texts:
Based upon a QGIS project I've created a simple webGIS using qgis2web plugin.
I've made some tests to improve the webGIS appereance using Google Chrome developer tool, but when I refresh or reboot the browser all edits went lost.

So I've edited CSS files with a text editor, then saved it and they became permanent. But when I open the index.html file with a text editor, it looks like doesn't contain page contents in the , while I expected to find the same as the one on Chrome console, but just javascripts.

How can I make permanent the Chrome development tool edits or edit directly the HTML file with a text editor?
I've got basic knowledge on HTML and CSS, but none in JavaScript.


texts: Can ArcPy's AsShape be used with GeoJSON FeatureCollection of Polygons with Labels?


texts: PostGIS ERROR: Relate Operation called with a LWGEOMCOLLECTION type

texts:
I have installed PostGIS on my machine and I am using shapefiles provided by the US Census Bureau for finding zips and counties of coordinates. I am just using the PostGIS provided method for loading data.

I am testing my query using this data: https://gist.github.com/erichurst/7882666

This is my query:

select
  main.latitude,
  main.longitude,
  c.cntyidfp,
  s.statefp,
  z.zcta5ce
from (
  select
    latitude::numeric(15,8) as latitude,
    longitude::numeric(15,8) as longitude,
    st_setsrid(st_makepoint(longitude::numeric(15,8), latitude::numeric(15,8)), 4269) as mypoint
  from temp.sample
) as main
left outer join tiger.county as c
  on st_intersects(c.the_geom, main.mypoint)
left outer join tiger.state as s
  on s.statefp = c.statefp
left outer join tiger.zcta5 as z
  on st_intersects(z.the_geom, main.mypoint)

Now, all of the test data and the production data I will be using are coordinates (single points), but I keep getting this error:

  ERROR:  Relate Operation called with a LWGEOMCOLLECTION type.  This is
  unsupported. HINT:  Change argument 1:
  'GEOMETRYCOLLECTION(LINESTRING(-71.255183 42.736524,-71.255147
  42.736554),POLY...'

It may be triggered by other points, this was only an example.

Why would a linestring be created? I am only testing points. I would like to be able to run the query on all of the data at once, not try each record individually. Is there something I can do to extract a point? Am I missing something? I am new to PostGIS but I have done a bit of research on this issue and others and I cannot seem to get around this.


texts: Pipe Network: Computing the inflows for systems with multiple reservoir (input source)

texts:
Let's say, for the sake of simplicity, we have the following Y Shape water pipe network:

S1 and S2 are the reservoirs ( input source), Sout is the water demand where water is drawn out for commercial use or whatever use. We know the height of the reservoirs at S1 and S2 relative to the ground and we also know the height of the Sout.

All S1, S2 and Sout are inm^3/s.

For the pipe, we know for each pipe their corresponding diameter d, length L, roughness coefficient C.

My question is, given the above information, is there anyway we can uniquely compute S1 and S2 for the Y shape network above? From what I know, we have one equation but two unknowns ( the conservation of source):

Sout=S1+S2

But this is not sufficient for us to uniquely determine S1 and S2 unless we have other equations here. I think the height of the nodes are also playing a role in the solution, but I just don't know how it can be done. Any ideas?

On a side note, if our water pipe network is a loop, we can use the hazen william formula and hardy cross method to compute the head loss. This means that whatever that is required to be known about the pipe is already known, as long as they are used in hazen william formula and hardy cross method.


texts: Automatically change the labelling style depending on the visibility of each label on a layer?

texts:
I would like to show as many  street name/label as possible on the print composer map, without overlapping, automatically change each label’s style, if it is necessary, depending on what is possible to display.

The default label is curved, 8p (Liberation Sans), but if the QGIS can't show a label (because of any reason), then the style of label of that street should change to a narrow font (Liberation Sans Narrow). If it is not possible at all to show a label in curve, then it should change to parallel label.

Summarize (for each label on the layer):

If it is possible to show: Curved Label, normal font
Else if it is possible to show: Curved Label, narrow font
Else if it is possible to show: Parallel Label, narrow font
Else: Parallel Label (just the first 3 character and a ‘.’ For example: “Main street” -> “Mai.”)

How is it possible in QGIS?


texts:
When I load a vector with a large number of features into QGIS, the rendering time can become extremely slow. But when I look at my system resources, they all have plenty of excess capacity.

What is the bottleneck that makes the rendering take so long? Is it a question of system resources that for some reason are not being allocated to QGIS?

If it matters I'm currently using an Intel i7 quad core / 16gb ram / NVIDIA gaming GPU / Windows 10.

What I'm trying to figure out is whether it may be possible to speed things up by upgrading hardware?


texts: How to extract max elevation values for each road links in a shapefile using Python libraries?

texts:
I want to extract values from a raster DEM file for each feature (polyline) in a shapefile. For each feature in a road shapefile, I want to extract the max and min elevation within a buffer (say 10 ft) of that line.

I am looking for a function similar to AddSurfaceInformation_3d from ESRI ArcGIS.

I attempted using GeoPandas and rasterio with no success. I am receiving only null values.

import geopandas as gpd
gdf = gpd.read_file(pth_to_road_shapefile)

from rasterstats import zonal_stats
with rasterio.open("path_to_test.tif") as src:
    affine = src.transform
    array = src.read(1)
    df_zonal_stats = pd.DataFrame(zonal_stats(gdf, array, affine=affine))

gdf2 = pd.concat([gdf, df_zonal_stats], axis=1)

Example DEM and vector shapefile here


texts: Comparing each geometry of shapefile with all geometries of another shapefile using ArcPy?

texts:
I want to compare each polygon from one shapefile with all the polygons from another shapefile.
When one polygon is geometrically identical to a polygon on the other shapefile, I want to copy the value of field "code" of the second shapefile to the field "code" of the first shapefile.
I wrote the following code:

fc1 = ur'D:\polygons1.shp'
fc2 = ur'D:\polygons2.shp'
fields = ['SHAPE@','code']
with arcpy.da.UpdateCursor(fc1,fields) as cur1:
    with arcpy.da.SearchCursor(fc2,fields) as cur2:
        for row1 in cur1:
            for row2 in cur2:
                row1[1] = row2[1] if(row1[0]==row2[0]) else ''
            cur1.updateRow(row1)

When I run the above code, the values from "code" are not copied plus there are no errors.
Any ideas?


texts:
How do I figure out the longitude and latitude coordinates necessary to create a circular 1 mile radius around one location? I don't mind how many coordinates that takes.

For instance: Latitude = 28.4789 Longitude = -81.4682

What mathematical theorem or formulas would you use to accomplish a task like this?

(I'm not using google maps, I just need to generate the necessary longitude and latitude to form the radius)


texts: What's the simplest way to style vector data in Geoserver?

texts:
There seem to be three main options:

Manually write the SLD files, based on the SLD cookbook
Create the SLD file in QGIS, then import
Use the add-on GeoExplorer to create custom styles online

Options 2 and 3 are attractive to me, as they avoid writing long SLD files. However, both have failed for me, meaning I still have boring grey polygons being served from my geoserver to a test map.

I also wonder if there other options I've not yet tried, such as styling on the client side, or creating styles with a different GIS whose .sld files work better with geoserver.

In visual terms, I want to go from this:

To this:

ASAP (as simply as possible!)


texts: Is there free desktop WFS client with BBOX support?

texts: To troubleshoot WFS server (geoserver) with big dataset (global OSM btw) I have tried number of free GIS desktops which work on Mac: QGIS (1.8, 1.9dev), uDig, gvSIG - none of them seems to support WFS request filter for currently visible area. Which makes them quite unusable for any bigger work datasets. QGIS used to have "Only request features overlapping the current view extend", which should be exactly what I'm looking for, but this seems to have disappeared.


texts: Is it possible with Torque to show only the months where there is information to show?

texts:
I have a date column that ranges from 1998 to 2013, but each year the specific dates are always from the same three months (March, April and May).
I would want to show only these months (where there is data) instead of every month of every year, because most of the time there is nothing to show. Is there any way to do this in CartoDB?


texts: Take closest Points in categories to a specific Location and do calculation

texts:
I have a points layer with 50+ points in it, and another point layer with only one point in it (a centroid). These points are categorised in 4 groups (labeled 0-3).
I used distance matrix to work out the distance from the points to the centroid, which are now in the points layer table.
I want to:
a) Take the value from the Observation column for the closest point for each of the 4 groups
b) I want to find the difference* in the value between group 0 and 2, and 1 and 4, so I end up with two values at the end.
*In excel this is calculated as ABS(a1-a2) for example.
An example table below:
Point   Observation Distance    Group
1   23  4   1
2   14  12  1
3   5   54  1
4   54  2   1
5   7   32  2
6   87  1   2
7   32  78  3
8   45  45  3
9   8   89  4
10  5   3   4

I've made some progress; I figured out that I can select the nearest point for each group by using Select By Expression:
Distance = minimum("Distance", group_by:="Groups", filter:= "Groups" = 1)
changing 1 for 1,2,3, 4 etc.
But then how to do the rest of the calculation?


texts: Adding a column to an attribute table leaves edit pencil icon grayed out in QGIS?

texts:
How do I add a column to an attribute table?

When I open up the attribute table I'm unable to click the edit pencil because it's grayed out?

Ultimately, I'm trying to transfer data in a column I have in Excel into this zoning shapefile.

I tried searching for a solution to this issue but wasn't able to find anything.


texts:
It is possible to set page size in ArcPy?

I'm looking for any tricky way to set the page size, using pure Python or ArcPy.

I'm using ArcGIS 10.1 for Desktop.


texts: Create sequential ID column type when loading GeoDataFrame to PostGIS

texts:
Is there a way to define GID column as sequential (auto-incremental) when loading GeoDataFrame (created from Shapefile) to PostGIS?
I use Geopandas and "to_sql" method:

geodataframe.to_sql(tablename, engine, if_exists='fail', index=False, dtype={'gid':INTEGER(), 'user_name':Text, 'sessions':Text, 'geom':Geometry('MULTILINESTRING', srid=4326)})

I could not find any valid SQLAlchemy data type. It seems that none exists. Other option I see is: do not load GID column but create it after uploading data to PostGIS.


texts: Saving only selected features after doing a selection by location in QGIS using Python Terminal?

texts:
I wrote the following script to save only the selected features after have a selection by location done. But I have got an error when the saving line is executed. I don't want to load these layers in the map canvas.

How could I do it right?

I'm using QGIS 2.18.17 and Linux Mint 19.

#
import processing

# Load vector layers
rio = QgsVectorLayer('/home/tiago/Documentos/PRIM_Mineracao/rio/mata/rio_doce.shp', 'rio_doce','ogr')
drenagem = QgsVectorLayer('/home/tiago/Documentos/PRIM_Mineracao/rio/mata/rio_mata_disol.shp','rio_mata_disol','ogr')

# Selection by Location
sel = processing.runalg("qgis:selectbylocation", drenagem, rio,['intersects'],0,0)

# Save selected features
crs = QgsCoordinateReferenceSystem("EPSG:4674")
nome1 = '/home/tiago/Documentos/PRIM_Mineracao/rio/mata/'
nome2 = 'rio_doce_P'
caminho = nome1+nome2+'.shp'
QgsVectorFileWriter.writeAsVectorFormat(sel, caminho, crs, "ESRI Shapefile", True)

That's the error message thar I've got:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: arguments did not match any overloaded call:
  writeAsVectorFormat(QgsVectorLayer, str, str, QgsCoordinateReferenceSystem, driverName: str = '', onlySelected: bool = False, errorMessage: str = '', datasourceOptions: Sequence[str] = [], layerOptions: Sequence[str] = [], skipAttributeCreation: bool = False, newFilename: str = '', symbologyExport: QgsVectorFileWriter.SymbologyExport = QgsVectorFileWriter.NoSymbology, symbologyScale: float = 1, filterExtent: QgsRectangle = None, overrideGeometryType: QgsWKBTypes.Type = QgsWKBTypes.Unknown, forceMulti: bool = False, includeZ: bool = False, attributes: Sequence[int] = [], fieldValueConverter: QgsVectorFileWriter.FieldValueConverter = nullptr): argument 1 has unexpected type 'dict'
  writeAsVectorFormat(QgsVectorLayer, str, str, QgsCoordinateTransform, driverName: str = '', onlySelected: bool = False, errorMessage: str = '', datasourceOptions: Sequence[str] = [], layerOptions: Sequence[str] = [], skipAttributeCreation: bool = False, newFilename: str = '', symbologyExport: QgsVectorFileWriter.SymbologyExport = QgsVectorFileWriter.NoSymbology, symbologyScale: float = 1, filterExtent: QgsRectangle = None, overrideGeometryType: QgsWKBTypes.Type = QgsWKBTypes.Unknown, forceMulti: bool = False, includeZ: bool = False, attributes: Sequence[int] = [], fieldValueConverter: QgsVectorFileWriter.FieldValueConverter = nullptr): argument 1 has unexpected type 'dict'
  writeAsVectorFormat(QgsVectorLayer, str, QgsVectorFileWriter.SaveVectorOptions, newFilename: str = '', errorMessage: str = ''): argument 1 has unexpected type 'dict'


texts: ArcMap Geoprocessing Tools Security Issue

texts: I have been having issues when running geoprocessing tools in ArcMap 10.6 where the tool crashes and refuses to popup fully. The issue can be resolved by changing the internet security setting from "medium-high" to "medium" but my companies IT/Security team does not approve. They are requesting the URL that the tools are trying to access but I have not been able to find any literature from ESRI about what exactly they are accessing. My IT team wants the URL so they can addd it to their trusted sites rather than change the internet security setting (which could open us up to other issues).


texts: Format of Date/Time Field in Datasets for use in CartoDB Torque CAT

texts:
I have a Date/Time field in my data set on cartodb that is formatted as follows: 2016-02-09T10:30:00Z  that is the value in the first row. The data type is set to 'date'. I want to use this field to create a Torque visualization where the slider on the map shows the date and time correctly. However, when I try this and select my Date/Time field, the slider that is added to the map (by default) is starting at 02/08/2016 1900

Is this a bug in Cartodb or the way I have my times formatted?

This is the link to my map: https://milam.cartodb.com/viz/f73a1f4a-c9e5-11e5-868e-0ef7f98ade21/public_map


texts:
Can someone explain why you need to delete the cursor objects and the iterator when you use for example:

arcpy.UpdateCursor(blah)

I've checked the esri docs and there doesn't seem to be an explanation there, but the object is shown to be deleted at the end of both example scripts.


texts:
I am trying to use updatelayer in three mxd's. Each mxd has 3 data frames. I am getting an error, can anyone help?

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

env.workspace = r"C:\Project"
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"C:\Project\\" + mxdname)
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    updateLayer = arcpy.mapping.ListLayers(mxd, "ways", df)[0]
    sourceLayer = arcpy.mapping.Layer(r"C:\Project\layers\ways.lyr")
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    for df in dfList:
        for lyr in arcpy.mapping.ListLayers(mxd, "", dfList):
            if lyr.name == "ways":
                arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, True)
                print 'UpdateLayer'
    mxd.save()
del mxd


texts: Using select by attributes then doing spatial join on selection in ArcPy?

texts:
I am using ArcMap 10.6.
How can I write the code to run spatial join on each match between Layer A and Layer B?
First step is to select 1st value in Layer A field (GUID) and select matching records in Layer B (field: GUID). Then do a spatial join from B to A, adding field (AVG), then moving to next value in Layer A and start over until end.
Trying to repurpose some code, but getting stuck...
import arcpy
arcpy.management.SelectLayerByAttribute(r"C:\evening.shp", "NEW_SELECTION", "GUID = ''", None)

arcpy.management.CalculateField(r"C:\drivetimeareas", "GUID", '"AVG", None)
print ("Done")```


texts: Problem with map tiles loading with Leaflet and Bootstrap

texts:
I want to create website with Leaflet map placed in bootstap tab. If I switch into tab where is Leaflet instantion tiles are not loading correct(I see the grey background). I guess there is a problem with size of a  because if resize the browsers window tiles load correctly. Any ideas?
Here is my basic concept.

https://jsfiddle.net/rjper9nz/6/


texts:
I am using PostGIS 2.0 and I want to use the raster functionality through mapserver. It seems as if the latest official GDAL release 1.9.1 has an unstable PostGIS raster driver, as pointed out in a previous question. After posting a ticket, was advised on github to update GDAL from trunk for this to work. I did that but I still get the same error. Thus, my question is:

Is it necessary to build PostGIS again from source after updating GDAL for that update to become effective in PostGIS?


texts: how to correct merge rasters with different projections?

texts:
how to merge rasters with different projections with correct transforn? For example  EPSG:32636 and EPSG:32637

if i do it with gdal_retile command, i have correct view for EPSG:32637 rasters and uncorrect for EPSG:32636


texts:
How can I create a 3D shapefile containing line features or point features?
I have to extract the elevations from a DEM and then export a shapefile containing 3D features.
I'm currently working with QGIS 2.8.4


texts: Execution Error in Python when calling Clip tool

texts:
I have written a simple code in python for clipping. The code is :

    #Step 1 To clip the Study area
    import arcpy
    import os
    import sys
    from arcpy import env
    from arcpy.geoprocessing import gp
    #Set Workspace
    arcpy.env.workspace = r"C:\test"
    # set local variables
    in_features = arcpy.GetParameterAsText(0)
    clip_features = arcpy.GetParameterAsText(1)
    out_feature_class = arcpy.GetParameterAsText(2)
    xy_tolerance = ""
    # Execute Clip
    arcpy.Clip_analysis(in_features, clip_features, out_feature_class,  xy_tolerance)

On execution I am getting this error message:

 Executing: Script C:\test\lithology.shp C:\test\study.shp 0.0001 C:\test\lithoclip.shp
 Start Time: Wed May 30 18:41:05 2018
 Running script Script...

Traceback (most recent call last):
File "D:\kiran\test.py", line 9, in <module>
import arcpy
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\__init__.py", line 21, in <module>
from arcpy.geoprocessing import gp
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\__init__.py", line 14, in <module>
from _base import *
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 578, in <module>
env = GPEnvironments(gp)
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 575, in GPEnvironments
return GPEnvironment(geoprocessor)
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 528, in __init__
import weakref
File "C:\Anaconda2\Lib\weakref.py", line 14, in <module>
from _weakref import (
ImportError: cannot import name _remove_dead_weakref
Failed to execute (Script).

Can anyone help me understand where I am making the mistake?


texts:
How I would add WELD satellite imagery from NASA into ArcMap?

I'd like to take the entire CONUS scan and place it into ArcMap.

I tested a smaller scan of the Northeast, however I got a plethora of separate images.


texts: How to plot geographical coordinates in mapinfo?

texts:
I want to get started with mapinfo and plot some coordinates i.e latitude and longitute
. How is it possibe?


texts: I have a single band RASTER of an area, in which I want to create 360 bands and then assign those bands values from a table.  Is there a process for doing this?


texts:
I have a map with graduated symbology, with 328 unique values. I would like to add a legend to my map (as in the picture), illustrating that light blue indicates low values and dark blue high values.

I have done this in ArcGIS before. Is there a way to do this in QGIS as well?


texts:
I can't import arcpy into 64 bit Python. It works just fine in the python install from arcGIS, though. When I try to import it elsewhere, I get this error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\__init__.py", line 21, in <module>
from arcpy.geoprocessing import gp
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing\__init__.py", line 14, in <module>
from _base import *
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 14, in <module>
import arcgisscripting
ImportError: DLL load failed: %1 is not a valid Win32 application.

I've been looking for a solution, but what I'm finding seems to be "too bad, doesn't work in 64-bit." Is this the case?

Edit: working with ArcGIS 10.2.2, Python 2.7.8

Edit: I believe the problem is that I only have the trial version of ArcMap, which apparently doesn't support the 64-bit background geoprocessing installation


texts:
I want to write a Python function to create a shapefile layer from clipping an input layer using a mask.

Does someone know the way how to do it?


texts: Removing layer from TOC when it has space in name using ArcPy?

texts:
I am trying to remove layers from the TOC using Python Toolbox. It works provided that the layer names have no spaces in their names. Tho whole script looks like:

class LayerprocessingRemoveUntickedLayers(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Remove unticked layers"
        self.description = "Removes unticked layers."
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        param0 = parameter("toclayers", "Unticked layers", "GPString", "Optional", multiValue=True)

        params = [param0]
        return params

    def updateParameters(self, parameters):
        unticked_set = set()
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if not lyr.visible:
                unticked_set.add(lyr.longName)
            del lyr

        parameters[0].filter.list = sorted(list(unticked_set))
        return

    def execute(self, parameters, messages):
        toclayers = parameters[0].valueAsText.split(";")  # ??? Here is my problem

        arcpy.AddMessage("\nTicked layers:")
        for t in toclayers:
            arcpy.AddMessage("... " + t)

        arcpy.AddMessage("\nRemoving chosen layers... ")
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if not lyr.visible and lyr.longName in toclayers:
                arcpy.AddMessage("... " + lyr.longName)
                arcpy.mapping.RemoveLayer(df, lyr)
        return

Step 1: updateParameters function reads layers and gives me layers´ names.

Step 2: I tick layers I want to remove (it does no matter whether they are broken or not). There is no group layer.

Step 3: the script removes only those chosen layers that have no spaces in their names (in this case PP_N_eroze_SUp and PP_N_eroze_SUt). Obviously because layers with spaces in their names get quotes to their values.

Mr. Che commented that layers with spaces in their names get quotes when passed from tool dialog to script executing. And now I remember that I asked a question with the same problem some time ago as Multivalue input from chosen field giving patchy input strings from ArcPy?

The problem was in arcpy.GetParameterAsText(0), the solution was changing it to arcpy.GetParameter(0) and splitting its strings with ";" as Hornbydd suggested a long time ago.

It worked when the tool was within an ordinary toolbox. Now I have this tool in the Python toolbox but changing parameters[0].valueAsText.split(";") to str(parameters[0].value).split(";") is not working. How to get values (this is to say split values from parameters[0]) correctly?


texts: How to make icons from Layers Panel in QGIS to be larger?

texts:
I need to change the size of layers icons  placed in Layers Panel of QGIS.
I have tried to find some solution inside Settings/Options but without success.


texts: Coordinates remain very high when converting projections

texts:
I have this GeoJSON file which uses the following projection: EPSG:28992. I want to convert this to EPSG:3857 but I don't have GIS software available. I have tried ogre, but when I first convert it to a shapefile and then back to GeoJSON using the aforementioned projections, it still shows coordinates with extremely high numbers (476628.5890141478,6820821.561640609 etc). The same applies when I try converting it using reproject.

In both cases when I try loading the converted files into geojson.io it fails -- it doesn't show on the map. In EPSG:3857, I'm expecting something like 4.30, 52.08 instead of coordinates that are in the hundreds of thousands.

Any idea what I am doing wrong or how to get the right values?


texts: Convert an imported vector FeatureCollection to raster image with correct longitudes and latitudes and a pixel scale of 30m

texts:
I have successfully imported a CSV file detailing pixel land surface temperature (LST) difference, longitude and latitude into GEE as a FeatureCollection. Please see: https://code.earthengine.google.com/77502e1432d8188f0d68f3fc37684adf
However, the data points projected onto the map are very small so I tried to increase the pixel size using:
var LST = LST.reproject({crs:'EPSG:4326', crsTransform: null, scale:30})
But now the pixel centres no longer match the initial longitudes and latitudes and the pixels don't actually seem to 30x30m, please see below:

How can I make the imported LST difference pixel centres match the imported longitudes and latitudes with 30x30m meter pixel sizes? Perhaps the data I am importing is in the incorrect format or lacking something.


texts: Creating features into a layer by using Model Builder

texts:
My friend created a ModelBuilder to place a geographic mean centered point between two selected road lines. Then it calculates the distance between road lines and types its value to a field. After that it calculates the distance by multiplying it with 2 and "boom!" here the distance between two road lines. Then she creates a buffer with this point to create a circle in another layer to show the distance value in it symbologicaly.

She wants to do this process continuously, i mean first she selects two road lines and run the model and she gets what she wants, then she wants to select another two road lines and get a new feature in this layer. But the problem is everytime she runs this code, it creates the feature (and layer) but when you run it again it deletes the previous feature. So everytime you run it, you get only 1 feature. Copying the created feature to another layer didn't work either.

Can you help us to make the model run and store the created features in a layer so she can have a layer full of road distances?

enter image description here


texts: Adding PostGIS query layer says Query class name is invalid?

texts:
I am having issues with adding postgis query layers in ArcMAP 10.2.2.3552, not doing anything special just adding a layer, this works in 10.2.1 so I don't know if its a bug or something

Also this error only happens when adding a query layer from a database over the network using ip address, the same database in the local install has no issues


texts:
I am trying to set the Index of a wms layer called "landsat1". It draws on top by default, and I want to move it down the order. I have tried

landsat1.setZIndex(1)

But I get an error "function not found". Im using ol 3.11.0. According to the api documentation, the function is experimental....but it doesn't seem to work for any version that has this function....other functions like landsat1.getProperties() works fine....why does this function not work? Or is there a better way to set the zindex om the wms layer?


texts: How to produce multiple maps with same layout for different layers' extents?

texts:
I have an area which is divided into several zones the borders of which I have defined as separate layers.  The layers are defined by sql queries on a postgis table of line segments joined with a table that says which segments are part of which zone boundaries.

What I want to do is produce a set of maps where each map is a zoom to layer extent for each zone and which has standard header, legend and scale (which will preferably vary with the zone but could be fixed across all maps).

My first attempt was to get a new composer for each map -- that worked but was rather tedious.  Next I tried producing one map and printing it and then deleting the map frame and going into qgis to zoom to the next layer extent and make the boundary visible then hopping back into the composer to load it in.  This worked more or less.  Map was fine, legend was OK because it was the same for all the maps and the zone name was easy to edit but the scale needed to be deleted and set up again for each map.

Next I looked at Atlas Generation.  The manual seems to assume that you know how it works and all you need is a reminder about the details (like so many software manuals ;)  It seems to assume that you have something that you are iterating over but I can't see how to make it work -- there is nothing to iterate over to get the layers.

Any suggestions welcome!!  :)

Russell


texts:
I am trying to do a simple buffer around a geojson I created with a model.  The buffer does not show up.  Here is my code.

    lyrStores2 = L.geoJSON.ajax("{% url 'stores' %}", {

    }).addTo(mymap2);
    storeBuffer = turf.buffer(lyrStores2.toGeoJSON(), 10, {units: 'miles'});
    lyrStoreBuffer = L.geoJSON(storeBuffer, {style:{
            color:'green',
        }
    }).addTo(mymap2);

The layer itself appears just fine.  I'm not sure if the geoJSON isn't appearing because I havn't created a url or a view for it like I had to for the shapefiles I created models for...

I have tried linking the  for the minified js and even tried this cdn

src="https://cdn.jsdelivr.net/npm/@turf/turf@5/turf.min.js"

Update:
So I saved the file as a geoJson using qgis.  I put the file into my static folder and everything worked fine.  This is what the code looks like for that:

    Storesjson = L.geoJSON.ajax("{% static 'mymap/data/storesgeojson.geojson' %}");
            Storesjson.addTo(mymap);

    Storesjson.on('data:loaded', function(){
            jsonStores = turf.buffer(Storesjson.toGeoJSON(), 10, {units: 'miles'});
            storesBuffer = L.geoJSON(jsonStores, {
                    style: {
                        color: 'red',
                    }
                }).addTo(mymap);
            });

The original file is loaded in like this

lyrStores = new L.GeoJSON.AJAX("{% url 'stores' %}", {

            });

My url for that looks like

path('Store_data/', Stores_datasets, name='stores'),

This is my view where I have it serialized:

def Stores_datasets(request):
stores = serialize('geojson', Stores.objects.all())
return HttpResponse(stores, content_type='json')


texts:
I had some problems with the processing alg reprojectlayer in QGIS3. Earlier in QGIS2.18 I used following command:

  processing.run('qgis:reprojectlayer', 'path_AND_filename_to_org_file','EGSG: 4326', 'path_to_new_file')

But it doesn't work in QGIS3


texts: Drawing line from variable number of points stored in two lists

texts:
I have two Python list latitude and longitude which stores the coordinates, and I want to draw a single line from the list coordinates.
I can draw the line if I know the number of latitude/longitude.
# suppose the line is having 2 longitudes and 2 latitudes.
QgsFeature().setGeometry(QgsGeometry.fromPolylineXY([QgsPointXY(78.3022187,23.7875547),QgsPointXY(98.3022187,33.7875547)]))

# suppose the line is having 3 longitudes and 3 latitudes.
QgsFeature().setGeometry(QgsGeometry.fromPolylineXY([QgsPointXY(78.3022187,23.7875547),QgsPointXY(98.3022187,33.7875547),QgsPointXY(148.3022187,133.7875547)]))

How can I draw a line in QGIS through Python, for the variable number of points?


texts: POST request on WFS is not working with URL but its working with IP address

texts:
I set up new workspace on GeoServer. Transactional WFS is enabled on that workspace. But I receives an error:
<?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://mywebpage.com/schemas/ows/1.0.0/owsExceptionReport.xsd">
    <ows:Exception exceptionCode="NoApplicableCode">
        <ows:ExceptionText>Transaction support is not enabled</ows:ExceptionText>
    </ows:Exception>
</ows:ExceptionReport>

I tested a POST request in postman and its working if I use IP adress instead of URL.
What can cause this problem to work with IP and not with URL?


texts:
I am trying to build raster summarizing features into a GeoDjango webmap,
and would like to use the django-raster and django-wms packages to do so. I believe that I have them set up correctly, but am hitting errors when I try to install the MapScript Python bindings for MapServer (necessary for
django-wms).

I believe there is a compatibility error with Python 3 and MapScript, but
am unsure of how to resolve them. I found this forum-post & GitHub issue , but unfortunately my software development knowledge is not solid enough to make much sense of them.

Has anyone compiled MapServer/MapScript with Python 3 recently, or know of any information per updating the software for compatibility with Python 3+?


texts: I have set up an operation dashboard to help keep track of field work. The field staff are using the field maps app to collect information and photos of point features. Is there anyway to get the photos to show up in one of the widgets? I can access them through pop ups in the map but viewing them larger opens a new tab and I would like to avoid that if possible.


texts:
The situation is this:

I need to remove the lines inside the blue polygon, I have tried with dissolve and delete holes processing tools but they didn't work.
I would like to know if there is a solution. File: https://drive.google.com/file/d/1itv2Gpr4fPN6TpeSKJDRCZEox0hnFP76/view?usp=sharing


texts: Baselayer disappears when I draw a layer

texts:
I have three layers, the base layer and two WMS layers. When I load the map, one of them always shows up even when I put visibility:false. And when it shows up, the base layer disappears. What can I do?

My code in OpenLayers is:

map = new OpenLayers.Map('map', {controls: [], displayProjection: new OpenLayers.Projection("EPSG:32632"), numZoomLevels: 31, tileSize: new OpenLayers.Size(256,256)});
wms1 = new OpenLayers.Layer.WMS('base', 'http://172.26.10.126/fcgi-bin/mapserv.exe?map=/ms4w/apps/guinea/guinea.map', {layers: 'base'},{isBaseLayer:true, visibility:true,transparent:true});
wms2 = new OpenLayers.Layer.WMS('guinea','http://172.26.10.126/fcgi-bin/mapserv.exe?map=/ms4w/apps/guinea/guinea.map', {layers: 'guinea'},{isBaseLayer:false, visibility:true,transparent:true});
wms3 = new OpenLayers.Layer.WMS('municipios','http://172.26.10.126/fcgi-bin/mapserv.exe?map=/ms4w/apps/guinea/guinea.map', {layers: 'municipios'},{isBaseLayer:false, visibility:false,transparent:true});
map.addLayers([wms1,wms2,wms3]);
map.zoomToExtent(new OpenLayers.Bounds(123360.47,-162161.43,767730.83,418254.41));
map.setCenter(new OpenLayers.LonLat(8.80899,1.83925), 7);
map.addControl( new OpenLayers.Control.LayerSwitcher() );

And my map file is:

LAYER
NAME "base"
STATUS ON
PROCESSING "LABEL_NO_CLIP=True"
CONNECTIONTYPE POSTGIS
CONNECTION "host=localhost user=postgres dbname=Guinea_Ecuatorial password=tph"
DATA "geom from gnq_africa_central using unique gid srid=32632"
METADATA
    WMS_NAME "base"
    WMS_SRS "EPSG:32632 EPSG:4326"
END
TYPE POLYGON
LABELITEM "na2_descri"
CLASS
    STYLE
        COLOR 185 122 87
        OUTLINECOLOR 0 0 0
    END
    LABEL
        TYPE TRUETYPE
        FONT "arial"
        SIZE 15
        COLOR 0 0 255
        ANTIALIAS TRUE
    END
END

LAYER
NAME "provincias"
STATUS ON
PROCESSING "LABEL_NO_CLIP=True"
CONNECTIONTYPE POSTGIS
CONNECTION "host=localhost user=postgres dbname=Guinea_Ecuatorial password=tph"
DATA "geom from gnq_provincia using unique gid srid=32632"
METADATA
    WMS_NAME "provincias"
    WMS_SRS "EPSG:32632 EPSG:4326"
END
TYPE POLYGON
LABELITEM "provincia"
CLASS
    STYLE
        COLOR 85 185 107
        OUTLINECOLOR 0 0 0
    END
    LABEL
        TYPE TRUETYPE
        FONT "arial"
        SIZE 15
        COLOR 255 0 128
        ANTIALIAS TRUE
        POSITION CC
    END
END
LAYER
NAME "municipios"
STATUS ON
PROCESSING "LABEL_NO_CLIP=True"
CONNECTIONTYPE POSTGIS
CONNECTION "host=localhost user=postgres dbname=Guinea_Ecuatorial password=tph"
DATA "geom from gnq_municipios using unique gid srid=32632"
METADATA
    WMS_NAME "municipios"
    WMS_SRS "EPSG:32632 EPSG:4326"
END
TYPE POLYGON
LABELITEM "distrito"
CLASS
    STYLE
        COLOR 85 185 107
        OUTLINECOLOR 0 0 0
    END
    LABEL
        TYPE TRUETYPE
        FONT "arial"
        SIZE 10
        ANTIALIAS ON
        COLOR 0 0 0
        POSITION CC
    END
END

The layer that always shows up is "distritos".


texts: Opengeo Suite - What to do after a change of IP address?

texts:
I have an Opengeo Suite installed on a server. This server had a static IP address with which I could access the whole suite.
Our internet provider changed the static IP and after this things are not working anymore. I can't load any map from online sources, neither Google layers or Bing layers are not showing up.
I also get the following error in Geoexplorer:
This page was unable to display a Google Maps element. Please contact the site administrator. If you are the administrator of this site, please check the JavaScript console or check the following page for troubleshooting: http://g.co/mapsJSApiErrors

I didn't use a Google API Key until now.

So my question would be, what do I need to do in the case of an IP change? Do I need to update the IP address somewhere in the config files or anywhere else ?


texts: Find shortest path between points avoiding shapefile?

texts:
I am an intermediate python developer with very limited GIS experience. I have a list of latitude/longitudes of a ship at sea. I need to join them into a polyline using the shortest distance between them to trace the path followed by the ship. However, the path cannot pass over land (for which I have a shapefile). How can I do this? I use Python 2.7.

I am presently not using any GIS software and hope to do the task using the python modules like matplotlib and scipy.


texts: Adjusting ArcMap Model to Create 1000 Iterations of Random Data Files with Spatial Join?

texts:
I have created a model that creates a random point data set for a certain area, conducts a spatial join of these points to polygons, calculates the frequency of a certain field in the joined table, and exports this frequency table to excel.

Ideally, I need to be able to run this random point model 1000 times so I can construct a distribution of the exported field frequencies to use as a null model for statistical testing (in another program).

How do I go about doing this without having to manually manipulate the model after every run?

I'm unsure of how to avoid 1) removing the spatial join in the model manually after every iteration, and 2) the exported excel being overwritten due to the output file having the same name as the previous iteration of the model.

If you're able to provide tools I can use within ModelBuilder, that would be preferred.

Current Model: input layer -> create random points (tool) -> output and polygon layer for join - > spatial join (tool) -> output -> frequency (tool) -> output -> table to excel (tool) -> excel file


texts: Method for estimating accuracy in spatial interpolation?

texts:
When using a spatial interpolation method (any method) without having a way to validate it, besides sampling some points out, is there a way to count as well with the mean distance between observations and spatial distribution.

I wonder if there is another way to give one quality/accuracy measurement value, instead of giving all three different measurements of the surface quality separately.


texts: Data not filtering with popup paging (ArcGIS JS API 4.15)

texts:
I have created and drawn a filtered layer on my map... Filter has worked fine with the drawn graphics; however I have just now noticed that after enabling pagination in the popup (unhiding it), the popup paginates through the data I had filtered? I am querying earthquake data..  In most of my popup scenarios the pagination is paginating through the data I have filtered.. for instance, I am drawing say 3 earthquakes after the filter on the map... map shows three earthquakes, start paginating in the popup and you are sifting through 80... pagination is not honoring the filter? Anyway I can handle this? What am I missing here?
     // original filters which are being honoroed everywhere except in popup paging...

        query.where = "eventTime >= CURRENT_TIMESTAMP - 30 AND updated >= CURRENT_TIMESTAMP - 30 AND grid_value > 2";
        query.outFields = ["*"];
        query.orderByFields = ["id"];
        query.returnDistinctValues = true;
        query.returnGeometry = false;

     // some other code

            let popupShake = {
                outFields: ["*"],
                "title": "Shake ID: {id}",
                "content": function(feature) {
                        let name_plc;
                        let dlv = document.createElement("div");
                        dlv.className = "popd";
                        dlv.innerHTML = `Magnitude: ${feature.graphic.attributes.mag}
 ID: ${feature.graphic.attributes.id}
 Event URL: View
 Last Update: ${new Date(JSON.parse(feature.graphic.attributes.updated)).toUTCString()}
Intensity: ${feature.graphic.attributes.grid_value}
 Event Time: ${new Date(JSON.parse(feature.graphic.attributes.eventTime)).toUTCString()}`;
                 return dlv;
                },
            };

            fl = new FeatureLayer({
                source: gras,
                objectIdField: "ObjectID",
                geometryType: "polygon",
                fields: [{
                    name: "ObjectID",
                    alias: "ObjectID",
                    type: "oid"
                  }, {
                    name: "id",
                    alias: "id",
                    type: "string"
                  },
                  {
                    name: "mag",
                    alias: "Magnitude",
                    type: "double"
                  },
                  {
                    name: "updated",
                    alias: "Last Update",
                    type: "date"
                  }, {
                    name: "eventTime",
                    alias: "eventTime",
                    type: "date"
                  }, {
                    name: "url",
                    alias: "url",
                    type: "string"
                  }, {
                    name: "grid_value",
                    alias: "grid_value",
                    type: "single"
                  }],
                  renderer: renderer,
                  popupEnabled: true,
                  popupTemplate: popupShake,
            });

            app.mapview.map.add(fl);
            fl.legendEnabled = false;

            app.mapview.whenLayerView(fl).then(function(layerView) {
                eqLayerView = layerView;
                eqLayerView.filter = {
                    where: "id = ''"
                };
            });
        });


texts: Crash/error in ArcPy using Create Random Points?

texts:
I'm trying to deal with an odd crash that occurs irregularly in arcGIS's arcpy, specifically under the create random points tool.

This occurs under ArcGIS 10.1 to 10.3.

During "CreateRandomPoints", the program attempts to create by polygon a number of points in the polygon, the number dictated by a field.  The problem I experience is that occasionally, upon attempting this line, python will immediately crash without giving an error code.

I accept that it's a known issue. I accept that there's nothing to be done about the crash itself, other than working around it.  the problem is, i can't seem to catch the error and work around it.

I have embedded the code chunk containing

    arcpy.CreateRandomPoints_management(wkspc,outshape,shape,,field)

inside of a try/except block.  I have set up a logic where hopefully the code would execute (and if necessary, fail) 5 times before exiting the program with an informative error.  However, it only runs once, and exits the program without stating the error.  (my try/except block isn't catching the error).

I have tried every error code I can think relevant, even a naked except block which to my knowledge should catch everything.

Is there some other means to catch the error? something i haven't tried yet?


texts: How to serve KML files from ArcGIS Server in AGS JSAPI 3.6

texts:
I have a javascript test page that adds a couple kml layers to a map and I can toggle them off and on.  I seem to get info popups for free with KML.  When I click on a point feature I get the data associated with it.

I set up an ARCGIS for Server because I want to start doing some geoprocessing.  I can add a map service to my javascript as such:

var layer = new esri.layers.ArcGISDynamicMapServiceLayer("http://mywebserver/arcgis/rest/services/Test/PA_Hospitals_2000/MapServer")

However, when I do this I seem to lose the popup ability that happened with KML.  So I thought I would create a KML link.  For some reason the javascript api will not load the KML link when it is on ARCGIS for Server.  If I browse to the URL it downloads the file but it will not display on the map in the javascript api.

This just seems wrong as you should be able to serve KML from ARCGIS for Server.  I would rather have all my files on the arcgis server, from a workflow standpoint as opposed to exporting to KML for point features to get popups and other layers as a service.

How do you either serve kml files from ARCGIS for Server or retain popup data on a published service?


texts:
I need to process some MapInfo data in QGIS 3.4, but I am unfamiliar with MapInfo. The data seems wrong when compared to a Google Satellite base map. I do not know the projection of the MapInfo layer, and this is the GDAL-Vector information output:
INFO: Open of 'C:/.../LAYER.MIF'
      using driver 'MapInfo File' successful.

Layer name: LAYER
Geometry: Unknown (any)
Feature Count: 9
Extent: (12.291721, 41.737178) - (12.641386, 42.011683)
Layer SRS WKT:
GEOGCRS["unnamed",
    DATUM["Arc 1950",
        ELLIPSOID["Clarke 1880 (modified for Arc 1950)",6378249.145326,293.4663076,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Name: String (50.0)

I have tried reprojecting the layer to the two projections mentioned (EPSG:9001 and EPSG:9122), but the data remains wrong. In an associated .WOR file, I found this command:
Set Map
  CoordSys Earth Projection 1, 5
  Center (12.35610885,41.95712217)
  Zoom 60.26835134 Units "mi"
  Preserve Zoom Display Zoom
  XY Units "m" Distance Units "mi" Area Units "sq mi"

I do not know whether this is relevant, and couldn't find the exact meaning of Earth Projection 1, 5 on MapInfo documentation. Is there a way to know this layer's projection?


texts: Breaking up existing watershed polygons by flow to individual tributaries

texts:
I have an existing hydroshed dataset (HydroBASIN lev12 basins from the HydroSHED dataset) that I would like to split up further, so that each tributary has its own basin delineated (see below).

I have tried using the watershed tool in ArcGIS to build these sub-basins, using the same DEM that was used to create the HydroSHED dataset.  This works fine, but the boundaries do not line up exactly with the HydroBASIN lev12 basins due to the DEM pixel size. I'm not sure how they smoothed these basin lines when they created the HydroBASIN dataset.

I really need to keep the same HydroBASIN lev12 basin outlines (gray), and then split these polygons further into their sub-basins (orange).  I've tried "burning" the lev12 basin lines into the DEM before running the watershed tool, clipping the DEM by the lev12 basin line before running watershed, using Extend Line or Align Feature to manually adjust the lines inside the lev12 basin outline...none of these allow me to keep the same lev12 basin outline.
Any ideas?


texts:
I recently downloaded the OpenZoomstack mbtiles maps from OS, following the query:
Launching high-detailed Ordnance Survey maps in QGIS 3x
but it looks like the QGIS can't open them properly.
The way how to load the mbtiles in QGIS I found here:
https://support.maptiler.com/i77-open-maps-in-qgis
but it looks like the file is completely invisible in the program.

I thought that problem might lay in the coordinates system. However, the provider set the following CRS:
EPSG:3857 - WGS 84 / Pseudo-Mercator - Projected
The same as I set in properties for my project.
The layer is visible in the layers box but not on the map.
Where can be the problem then?


texts: Adding XYZ tiles to images with layer in QGIS

texts:
I have made a python script that loads several KML files and I'm gonna render them with a QPainter on a QImage. I set the extent to the QgsMapSettings.fullExtent(), when all my layers are part of the QgsMapSettings and added via QgsMapSettings.setLayers(layers). I only extent the scale of the view by 1.1 or 1.2 to just have a better view of the data with
rect = QgsRectangle(ms.fullExtent())
rect.scale(1.1)
ms.setExtent(rect)

This works all great and I gonna save the image with the QImage.save function. The only problem I have is, that I have a XYZ Tiles layer, based on OpenStreetMap and I'm not able to add it to my Image.
I tried picking the layer with QgsProject.instance().mapLayersByName("OpenStreetMap")and adding it to my layers list and finally add it to QgsMapSettings with the QgsMapSettings.setLayers(layers) method, but all that happens is that it shows the geo-coordinates as points of the KML files, but not the XYZ layer.
When I try to add the layer after my KML layers, just with QgsMapSettings.setLayers(osm_layer), then my image only shows the lowest zoom level, without any data on it, which makes sense, because all the layers, which I have added before, are gone.
So what exactly I'm doing wrong?


texts:
I am very new to PostGIS/PostgreSQL. In fact, I downloaded it because of the answers I got to this question: Distance between centroid and farthest point of polygon

I installed both PostgreSQL and PostGIS. But I am not able to import a shapefile to the database? It gives an error, "Shapefile import failed". What am I doing wrong?


texts: Run Arcmap commands from another application

texts:
I have a windows forms application (Z application), is it possible to start Arcmap, and start the Editor or any other command I desire from the application (Z appliction) ?

So far I was able to initialize and bind the license, so now I am able to access feature classes using Arcobjects without starting arcmap, but I want to start Arcmap, add featurelayers, start editing, and then save those edits. The only missing part is the part where I am able to run arcmap and access the application to get the document.

here is the flow

        RuntimeManager.BindLicense(ProductCode.EngineOrDesktop);
        IAoInitialize s = new AoInitializeClass();
        s.Initialize(esriLicenseProtCode.esriLicenseProductCodeAdvanced);

next step would be to start arcmap

        `Process.Start("arcmapDirectory")

Arcmap starts, then this is where I am stuck, I want to access the arcmap application, so I can get  mxdocument, and commands.


texts: Ranking occurrences in field using ArcPy?

texts:
I have a field C_CODE and I would like to count occurences of values in that field that will start with zero in front of the one digit numbers:Ex: 01,02 and normally for two digit numbers Ex:..09,10,11,etc

Wanted output:

C_CODE RANK  C_CODE_FINAL                        #Info#
17      01       1701      # RANK counts how many times each value in C_CODE exists
15      01       1501      # C_CODE_FINAL is a concatenation of the other fields
23      01       2301
29      01       2901
15      02       1502
19      01       1901
17      02       1702
15      03       1503
37      01       3701

and then concat the fields to C_CODE_FINAL .

Tried:

fc = path to shp
field = "C_CODE"

#Create dictionary to store unique values
count_values = {}
arcpy.AddField_management (fc, 'rank', 'LONG')
with arcpy.da.SearchCursor (fc, field) as cursor:
    for row in cursor:
        if not row[0] in count_values.keys():
            count_values[row[0]] = 1
        else:
            count_values[row[0]] += 1

It returns a rank column filled with zeros although there are many values that exist multiple times. Also they have to start with 01 for values that appear for the first time.


texts: R/ggplot2 adding an additional line/polygon to a choropleth map

texts:
I have a choropleth map that represents the number of states being affected by a storm event. This is shown in the image below. I want to lay a transparent line or polygon over these states to represent the storm area. How do I do this?

Here is the code I'm using for ggplot:

ggplot() + geom_polygon(data = final.plot, aes(x = long, y = lat, group = group, fill = factor(T3)), color = "grey", size = 0.3) +
          coord_map("conic", lat0 = 30, xlim = c(-121, -72), ylim = c(25, 50)) + theme(axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank(), legend.position = "bottom") +
          geom_point(data=s1.lat.long, aes(x=long, y=lat), color="black", size=1)

Here is a link to the whole R code and shape files.

EDIT: I would like it to look like this mock image below, whereby a band spreads all the way across the US, representing geomagnetic storm activity. This band can then be made fairly transparent.

Really, this is a crude way of reproducing this image:


texts: Maintaining extra data for .laz for use in R and Python?

texts:
I am working on releasing a benchmark dataset for trees in RGB + Lidar and I'd like the format to match the most common use cases and packages. I had been intending on proving a tree index using the user_data column. Alot of preprocessing happens in python and i'm writing from laspy. I realized that user_data is clamped 0-255 (there are thousands of trees), meaning that trees do not have a unique index. I'm starting to write an extra dimension function

https://stackoverflow.com/questions/50815580/appending-an-index-to-laspy-file-las

but looking here, it looks like lidR is going to ignore those columns

https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAS-class.html

and i'll need users to do lasadddata anyways. Perhaps it would be best to just provide a csv lookup table based on x,y,z? With millions of points, will this be the best strategy? I'm worried about precision rounding among platforms. Thoughts welcome on how future users will best interact with extra dim data.

EDIT: The question was deemed too broad. So here is a literal example.

Here is a python pandas dataframe (taken from a laspy-like object) with the x,y,z coordinates of a point and an associated label. What is the best way to write this label information and load it in the R lidR package.

pc.head()
              x            y       z   label
272  315547.689  4094399.467  12.288   562.0
287  315541.905  4094400.774   4.151  2233.0
289  315541.279  4094400.842   4.166  2233.0
291  315541.725  4094400.506   8.189  2233.0
292  315540.776  4094400.865   5.199  2233.0

One option is to write to csv and the use lasadddata function add the extra attribute "label" in R separately. However, for a reasonably large point cloud, this would create a 100MB file, which would then be read into R and then joined to ensure the vector order lined up

If you directly write the laspy point cloud, and stick the label data in the user_data column, values above 255 will be clamped by laspy. See link above.

ANSWER:

Because this was flagged, I can't answer my own question. As @JRR noted, the key aspect is to be mindful of the datatype. lidR will read the extra attribute. I hope this helps someone in the future.

So in python

import laspy

def write_label(point_cloud, path):

    #Create laspy object
    outFile1 = laspy.file.File(path, mode = "w",header = inFile.header)

    #First define the new types

    outFile1.define_new_dimension(
        name="label",
        data_type=5, #Data types may matter
        description = "Integer Tree Label"
     )

    # copy fields from laspy object
    for dimension in inFile.point_format:
        dat = inFile.reader.get_dimension(dimension.name)
        outFile1.writer.set_dimension(dimension.name, dat)

    outFile1.label = <label column from pandas here>
    outFile1.close()

lidR sees a new integer column "label"

> colnames(a@data)
 [1] "X"                          "Y"                          "Z"                          "gpstime"
 [5] "Intensity"                  "ReturnNumber"               "NumberOfReturns"            "ScanDirectionFlag"
 [9] "EdgeOfFlightline"           "Classification"             "Synthetic_flag"             "Keypoint_flag"
[13] "Withheld_flag"              "ScanAngle"                  "UserData"                   "PointSourceID"
[17] "R"                          "G"                          "B"                          "reversible index (lastile)"
[21] "label"


texts: About Spatial ETL and how much is it used in the GIS industry

texts:
Because I've been noticing that some internships are requesting some know-how of this tool.

I would like to ask how is it used in the industry?
What do you use it for at your job?
Where could you learn it?

as far as I know my Geomatics program doesn't offer courses in any related SpatialETL software.


texts: QGIS API to determine number of walls an AP intersects

texts:
We are trying to identify QGIS API which can determine where all does a ray of signal intersects in a given room. See image for reference.

What we are looking for as result is "the lat lng of all the points on the wall where the ray intersects the wall".
Any ideas / pointers will be great help.


texts: Leaflet styling geoJSON object radius by its feature property

texts:
I am trying to style geoJSON objects by its feature properties (Capacity). I tried 2 ways:

1:

function getWidth(Capacity) {
return Capacity >= 40 ? 12 :
       Capacity >= 30 ? 8 :
       Capacity >= 20 ? 6 :
       Capacity >= 10 ? 3 :
       0;
       }

function getSize(h) {
return h >= 40 ? 12 :
       h >= 30 ? 8 :
       h >= 20 ? 6 :
       h >= 10 ? 3 :
       0;
       }

  function mystyle(feature) {
  if (currentZoom <= 12) {
    return {
      weight: 1,
      opacity: 0.5,
    };
  }
  else {
    return {
      weight: getWidth(feature.properties.Capacity),
      radius: getSize(feature.properties.Capacity)
    };
  }
}

var johnlewis = L.geoJSON(jl, {
 style: mystyle,
 pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng);
},
  }).addTo(map);

and as a result I got something like this:

where only weight changed, although the zoom level 12 set in function (map.getZoom() doesn't work at all.

In my 2nd example:

  function style2(feature) {
return {
    fillColor: setColor(feature.properties.Type),
    fillOpacity: 0.5,
    weight: 2,
    opacity: 1,
    color: '#ffffff',
    dashArray: '3'
};
}

var johnlewis = L.geoJSON(jl, {
style: function(feature) {
                    if (feature.properties.Capacity >= 5) {
                        return {
                            weight: 1,
                            radius: 4
                        };
                    }
                    if (feature.properties.Capacity >= 10) {
                        return {
                            weight: 2,
                            radius: 6
                        };
                    }
                    if (feature.properties.Capacity >= 20) {
                        return {
                            weight: 3,
                            radius: 8
                        };
                    }
                    else style2
                },

    pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng);
    },
  }).addTo(map);

It looks slightly better, because the records has been distinguished both by weight and radius, however there is no default situation shown, that was descriubed in the else value. Although only 1st if responded. Another ones didn't contribute into this view, hence I have only 2 types of radius instead of 4.

Does anyone know how to solve this issue?

I want to have a full styled geoJSON by properties, but now I am suffering for a lack of default style.

I saw a few similar examples, but basically they doesn't match my situation.

EDIT:

Recently I tried stuff like this:

 function mystyle(feature) {
 var rad;
 if (feature.properties.Capacity = 3) {
   rad = 4;
    } else if (feature.properties.Capacity = 5) {
    rad = 6;
    }
    else if (feature.properties.Capacity = 18) {
    rad = 8;
    } else {
    rad = 1;
    }
    return {
    fillColor: '#ff4d4d',
    fillOpacity: 0.5,
    weight: 0,
    opacity: 1,
    radius: rad,
    color: '#ffffff',
    dashArray: '3'
   };
}

 var johnlewis = L.geoJSON(jl, {
 style:mystyle,
    pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng);
    },
   }).addTo(map);

but instead of the radius I have got changes with the border weight (weight).
See the result on the map below:

So it looks like the radius changes aren't supported in Leaflet are they?


texts: How to pass the Basic Auth while calling GETMAP request on GeoServer using OpenLayers

texts:
How can I pass the GeoServer credentials while adding an layer using WMS. When I am doing GETMAP request I am getting the error in Chrome

[Deprecation] Subresource requests whose URLs contain embedded credentials (e.g. https://user:pass@host/) are blocked. See https://www.chromestatus.com/feature/5669008342777856 for more details.

And When I tried in Firefox I am getting ERROR Code 403. So that is very clear that I can't pass the credentials in the URL while calling GETMAP Request. My code is:
function loadCustomURL() {
            var url = "https://user:pass@securewms.url.com:8443/geoserver/ows?service=wms&request=getmap&version=1.1.0&layers=some_layer&srs=epsg:3857&bbox=some_bbox&format=image/png&width=500&height=600"
            var source = new ol.source.TileWMS({
                url:url,
                serverType: 'geoserver',
                params: {
                    'LAYERS': '2-0',
                    'TILED': true
                }
            })

            source.on("tileloadstart",  function () {
                $('.message').html("Drawing..");
            })
            source.on("tileloadend",  function () {
                $('.message').html("");
            })
            source.on("tileloaderror",  function () {
                $('.message').html("Error Please try again");
            })
            var layer = new ol.layer.Tile({
                source:  source
            });
            map.addLayer(layer);
        }

For GETCapablities when I had called I am getting its blocked because of CORS Policy
My code:
function WMSLayerCall() {
            var parser = new ol.format.WMSCapabilities();
            fetch('https://user:pass@securewms.url.com:8443/geoserver/ows?service=wms&request=GetCapabilities').then(function(response) {
                return response.text();
            }).then(function(text) {
                var result = parser.read(text);
                console.log(result);
                var extent = result.Capability.Layer.Layer.find(l => l.Name === '9').EX_GeographicBoundingBox;
                console.log(extent);
                //map.getView().fit(extent, map.getSize());
                var water_layer = new ol.layer.Tile({
                    source: new ol.source.TileWMS({
                    url:"url",
                    serverType: 'mapserver',
                    srs: 'EPSG:4326',
                    bbox: extent,
                    format : 'image/png',
                    params: {
                        'LAYERS': "9",
                        'TILED': true
                    }
                })
            });
            map.addLayer(water_layer);
            });
        }

So how can I pas the credentials for a GETMAP request?


texts: ArcGIS Portal: Webadaptor cannot connect to Portal?

texts:
I'm trying to setup an ArcGIS for Portal 10.2.2 installation on Windows Server 2012, which are both running on a VM (Hyper-V). I have both installed Portal for ArcGIS as well as the Webadaptor and all i need to do is to connect them. The Portal Service is running and the website is accessible under http://localhost:7080/arcgis/home. The Web Adaptor is running with tomcat 7 and SSL. I have set up a self-signed SSL certificate. The connector in my server.xml looks like this:

<Connector port="80" protocol="HTTP/1.1"
               connectionTimeout="20000"
               redirectPort="443" />

<Connector port="443" protocol="org.apache.coyote.http11.Http11Protocol"
           maxThreads="150" SSLEnabled="true" scheme="https" secure="true"
           clientAuth="false" sslProtocol="TLS"
           keystoreFile="C:\Program Files\Tomcat 7.0\webadaptor.ks"
           keystorePass="passwort"/>

When i open the web interface of the adapter via http://localhost/arcgis/webadaptor a dialog shows up where i can choose between ArcGIS Server and Portal. When i choose "Portal" and click on continue i get an error saying that the url of the web adaptor configuration is only accessible from the machine that hosts the web adaptor. I don't understand that because it is the same machine and i'm accessing it via localhost.

Can someone explain me whats wrong, please?


texts:
I have a points shapefile, and a file that tells me which of the points have a line between them. I am using arcpy to create a new shapefile to hold the lines and then create the lines based on the specifications in the file. For creating the lines, I am reading the coords of the points from the point shapefile itself.

out_path = arcpy.env.workspace
out_name = arcpy.ValidateTableName("the_lines")
geometry_type = "POLYLINE"
has_m = "DISABLED"
has_z = "DISABLED"
spatial_reference = arcpy.Describe(points_shp).spatialReference # Use Describe to get a SpatialReference object
# Execute CreateFeatureclass
arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type,spatial_reference = spatial_reference)

array = arcpy.Array()
point = arcpy.Point()
tie_shp_file = arcpy.env.workspace + "/" + out_name
insert_cursor = arcpy.InsertCursor(tie_shp_file)
feat = insert_cursor.newRow()

for line in all_lines:   #all_lines is of the form array([[0, 3],[0, 7],[1, 3]) indicating line [from, to]
    #from
    point.x=coords_pts[str(pt[line[0]])][0]
    point.y=coords_pts[str(pt[line[0]])][1]
    print "FROM: " + str(pt[line[0]]) + " COORDS_X: " + str(coords_pts[str(pt[line[0]])][0]) + " COORDS_Y: " + str(coords_pts[str(pt[line[0]])][1])
    array.add(point)
    #to
    point.x=coords_pts[str(pt[line[1]])][0]
    point.y=coords_pts[str(pt[line[1]])][1]
    print "TO: " + str(pt[line[1]]) + " COORDS_X: " + str(coords_pts[str(pt[line[1]])][0]) + " COORDS_Y: " + str(coords_pts[str(pt[line[1]])][1])
    array.add(point)
    print "ARRAY: " + str(array)
    polyline=arcpy.Polyline(array)# Create a Polyline object based on the array of points
    array.removeAll() #Clear the array
    feat.shape = polyline
    insert_cursor.insertRow(feat)

coords_points is a dictionary that has the form:

 {'10': (8594242.514466941, 3330642.8703512326)...}
i.e. {'Oid_str': (x,y)}

The code is creating the required number of entries in the shapefile but no lines are being shown and the lengths of the line is shown as 0.

Where am I going wrong?


texts: How to batch export rasters as GeoTIFF with new CRS?

texts:
I need some help concerning the exporting of a large batch of raster-files (mostly SID-types) as GeoTIFFs in a reference system (WGS 84) that is different from the one in which the original files were georeferenced (Belgian Lambert 1972).

I didn't want to export the files one at a time, as I have over 2000 different files to go through. For that reason, I attempted to use my little knowledge on Python to come up with a script in QGIS that could automatically export a large batch of raster-files (see below).

This code worked fine at first, but afterwards I discovered in ArcGIS that the extent of my GeoTIFFs was still given in meters instead of degrees, even though the reference system in the properties clearly states WGS84 (see image). Can any of you tell me how I should adapt my code so that the extent is converted as well?

#Directory where the exports have to end up in
myDir = '//fsarchive/INBO_historische_kaarten/_EXPORT_/'

#Reference System that I want the exported tiffs to have
trs = QgsCoordinateReferenceSystem()
trs.createFromId(4326)

layers = iface.legendInterface().layers()
pipe = QgsRasterPipe()

#for every layer in the layer interface
for index, layer in enumerate(layers):
    name = layer.name()
    extent = layer.extent()
    width, height = layer.width(), layer.height()
    renderer = layer.renderer()
    provider = layer.dataProvider()
    pipe.set(provider.clone())
    pipe.set(renderer.clone())

    print (index + 1)   #these prints are just to follow the progress when the script is running
    print ("of")
    print (len(layers))
    print (name)

    opts = ["COMPRESS=ZIP"] #compress the GeoTIFF with the ZIP-method

    file_writer = QgsRasterFileWriter (myDir + name + ".tif")
    file_writer.setCreateOptions(opts)
    file_writer.writeRaster (pipe, width, height, extent, trs)


texts: Using Pgrouting dijkstra algorithm in PostGIS?

texts:
I own two tables one with streets and one points, i use the table with the roads to calculate routes with dijkstra algorithm. But I want when someone add one point in the table Point and this point is part of road from streets table ,then how i make the dijkstra algorithm not to take into consideration this road wich contain this point, when it calculate a route.

If i want to go from A --> B is give me the read way the dikjstra . i want if i put a point like the point of image the dijkstra change the way and give me the orange way. Can i make this?


texts: Import Cartesian coordinates from Excel and count its area

texts:
I'm new and struggling with calculating xy area from Excel data ( or it can be any type of data with xy coordinates). Moreover according to task it must use all existing coordinates in Excel file.

My code

import openpyxl

book = openpyxl.load_workbook('coordinates.xlsx')
sheet = book.active

for row_i in range(2, sheet.max_row + 1):
    x = sheet.cell(row=row_i, column=1).value
    y = sheet.cell(row=row_i, column=2).value
polygon = [[x,y]]
print(polygon)
area = 0.0
n = len(polygon)
for i in range(n):
    i1 = (i+1)%n
    area += polygon[i][0]*polygon[i1][1] - polygon[i1][0]*polygon[i][1]
area *= 0.5
print('area = ', area)

I don't know how to make that polygon = [[x,y]] would be like polygon [[608....,5105...],[6085.....,517...]], my code only picks last row data, I need that it pick all data and use it to calculate. Or if you have better idea, please share it.


texts: How can I create a layer for each unique field item in a shapefile using Arcpy

texts:
#Introduction and imports
import arcpy, string
from arcpy import env

#User Inputs: layer, field name, output workspace (optional)

strTargetLayer = arcpy.GetParameterAsText(0)
strFieldName = arcpy.GetParameterAsText(1)
strSuffix = arcpy.GetParameterAsText(2)
strOutputWorkspace = arcpy.GetParameterAsText(3)

TargetLayer = arcpy.mp.Layer(strTargetLayer)
TargetLayerDQ = TargetLayer.definitionQuery
Field = arcpy.ListFields(strTargetLayer, strFieldName)[0]
env.workspace = strOutputWorkspace
env.overwriteOutput = True
mxd = arcpy.mp.ArcGISProject("CURRENT")
df = df.listMaps("Drill Plan Map-Main")[0]
cursor = arcpy.SearchCursor(strTargetLayer)
lstValues = []

#Creating set(lstValues)
for row in cursor:
    lstValues.append(row.getValue(strFieldName))
lstValueSet = sorted(set(lstValues), reverse = True)
del lstValues

#Creating layers in map
arcpy.AddMessage("Creating layers in map...")
LyrList = []

for i in lstValueSet:
    try:
        if TargetLayerDQ !="":
            if i is None:
                ilyr = arcpy.MakeFeatureLayer_management(TargetLayer, "NULL" + strSuffix).getOutput(0)
                ilyr.definitionQuery = TargetLayerDQ + ' and ' + "\""+Field.name+"\"" +' is NULL'
            elif Field.type == 'String':
                ilyr = arcpy.MakeFeatureLayer_management(TargetLayer, i + strSuffix).getOutput(0)
                ilyr.definitionQuery = TargetLayerDQ + ' and ' + "\""+Field.name+"\"" +' = ' + "\'"+ i +"\'"
            else:
                ilyr = arcpy.MakeFeatureLayer_management(TargetLayer, str(i) + strSuffix).getOutput(0)
                ilyr.definitionQuery = TargetLayerDQ + ' and ' + "\""+Field.name+"\"" +' = ' + str(i)
            arcpy.mapping.InsertLayer(df,TargetLayer,ilyr,"AFTER")
        else:
            if i is None:
                ilyr = arcpy.MakeFeatureLayer_management(TargetLayer, "NULL" + strSuffix).getOutput(0)
                ilyr.definitionQuery = "\""+Field.name+"\"" +' is NULL'
            elif Field.type == 'String':
                ilyr = arcpy.MakeFeatureLayer_management(TargetLayer, i + strSuffix).getOutput(0)
                ilyr.definitionQuery = "\""+Field.name+"\"" +' = ' + "\'"+ i +"\'"
            else:
                ilyr = arcpy.MakeFeatureLayer_management(TargetLayer, str(i) + strSuffix).getOutput(0)
                ilyr.definitionQuery = "\""+Field.name+"\"" +' = ' + str(i)
            arcpy.mapping.InsertLayer(df,TargetLayer,ilyr,"AFTER")
        LyrList.append(ilyr.name)
    except:
        if arcpy.GetMessages(2) != "":
            arcpy.AddMessage("**FAILED TO EXECUTE FOR " + str(i) + "**")
            arcpy.AddMessage(arcpy.GetMessages(2))
        pass
#Saving layers to folder
if strOutputWorkspace != "":
    arcpy.AddMessage("Saving layers in target location...")
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name in LyrList:
            arcpy.SaveToLayerFile_management(lyr.name,lyr.name,"ABSOLUTE")


texts: Grouping polygons and other geometrical shapes

texts:
My situation:

I have administrative areas marked with polygons
Many areas have many, many, many islands
I need main area and all it's islands to change color when any of the islands or main is hovered
I currently need to make if statement for main area + each island in mouseover and mouseout event which would be a lot of if statements in total

Is there a solution or workaround for this?

Could I group them somehow and use that group or something to easily change options or do something else for the whole group?


texts: What's a good library for displaying many small maps on a web page

texts:
I would like to create a visualization on the web similar to the following:
http://namemapper.babynamewizard.com/namemapper/

Are there any good libraries for doing this?  Is OpenLayers light-weight enough to handle many maps at once?


texts: QGIS does not launch after plugin and python loading

texts:
I have a big problem with lauching QGIS. I have been using QGIS in the past, but now I can't launch it. I tried re-install, but it doesn't help. I tried install different versions of QGIS (1.7, 1.8, 2), but in every version process is finished after loading plugins and python...

When I want to launch QGIS via command line... QGIS is open, and after loading plugins/python nothing happens... Windows is closed, processed is finished.

In install post reports there aren't any errors.

Do you have any idea why?


texts: I'm using QGIS 2.14.7 and I want to change the  QGIS splash screen into something more suitable for my project. I noticed that in older versions there was a .png that could be changed, but the path (images/splash/splash.xcf.bz) doesn't exist in version 2.14.7 anymore. Does any of you know how to change it now?


texts:
Is anyone aware of a dataset for Canada that is similar to http://www.broadbandmap.gov/?
I contacted Connecting Canadians (https://www.ic.gc.ca/eic/site/028.nsf/eng/h_00587.html) to confirm that they currently don't have any published map data besides the images on their site. They indicated that the data might be available after all of the contribution agreements have been signed by the various providers, but this process's timeline goes until 2019.


texts: I am looking for a reliable hosted Overpass server where I can pay for N requests per sec/hour/day. I looked at https://switch2osm.org/providers/, and https://wiki.openstreetmap.org/wiki/List_of_OSM-based_services , and googled, but did not find any. I found several for geocoding (Nominatim), but nobody offers overpass API. Public servers have low limits for use and not very reliable.


texts:
I have a MS excel table with several thousand numbers that I need to select out of a single shapefile field with 20,000 numbers.
How do I build a sql statement in select by attributes to do this?


texts: PostgreSQL Query which returns geom column as WKT set in function

texts:
Given: We are accessing a PostgreSQL DB from a VB NET application. When the query completes the DB rows are delivered in a VB NET DataTable. The Postgre DB geom column is in the format of System.Bytes[] (a byte array).  I can convert the Byte array to the WKB format but converting the WBK to WKT format appears more complex. Hence, I was thinking I could return the geom col results as WKT.

My issue with this is I cannot use a SELECT * in the query to the DB.  It seems I will need to parameterize ALL the columns in the SELECT command since I want the geom column in a DIFFERENT format.

Obvious Query:

SELECT * FROM temporary.summmarytable WHERE conditions = 'xxxxx'

has to become:

SELECT col1, col2 , col3. col4....ST_AsText(geom) WHERE condition = 'xxxxx'

I'm not a fan of long paramterized column lists as shown above.  Is there a way to reset the geom output within a PostgreSQL function.