texts:
I am trying to rasterize a polygon in R using teh rasterize function but the output is not correct. I am comparing the result with the output produced by the rasterize tool available in QGIS

May data are:

> Fields2
class       : SpatialPolygonsDataFrame
features    : 1293
extent      : 219807.1, 230535.3, 4097488, 4104080  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 4
names       :  ID,       area, id_2,      mean
min values  :   1, 100044.627,    1, 1001.9653
max values  : 999,   9990.943,    1,  998.5807

> test
class       : RasterLayer
dimensions  : 660, 1074, 708840  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 219800, 230540, 4097480, 4104080  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : fortesting
values      : 126.7778, 1965.884  (min, max)

I am using the following command poly_to_raster<-rasterize(x=Fields2, y=test, field="mean")

The results looks like this:

The white pixels do not have any value and are NA. Confirmed by looking at the pixel info of this raster in QGIS

I thouhgh there is someting wrong wity the shapefile I am using as reference in rasterize, but when I do the same process in QGIS, I get the right result:

--- EDIT ---

The mean attribute comes from an extract function

meancalc<-extract(S_images_ts[[1]], Fields, fun=mean)
Fields@data<-data.frame(Fields@data, mean=meancalc)

For confirmation

> str(Fields2$mean)
 num [1:1293] 275 1350 442 430 1364 ...


texts: Managing layers from a specific group with PyQGIS

texts:
I have some layers loaded in QGIS, in various groups. I want, with PyQGIS, do some stuff on layers from one of theses groups, like adding geometry type to the layer name.
This way, I can do it on all layers :
layers = QgsProject.instance().mapLayers().values()
for layer in layers:
   layer.setName(layer.name()+"_"+QgsWkbTypes.displayString(layer.wkbType()))

This way, I can access layer names from one specific group :
root = QgsProject.instance().layerTreeRoot()
my_group = root.findGroup('My group')
for layer in my_group.children():
   print(layer.name())
   if isinstance(layer, QgsLayerTreeLayer):
      layer.setName(layer.name()+"_"+QgsWkbTypes.displayString(layer.wkbType()))

The "print" of layer name works, but I can't change name, either get the geometry type. I get this error :
AttributeError: 'QgsLayerTreeLayer' object has no attribute 'wkbType'

Does someone know what is the correct way to do that?


texts: Automate Splitting Polygons Into Sections

texts:
I need to split many cemetery sections into 8 (mostly) equal pieces each. Each section is for the most part square, but there are a few odd shapes too. And the squares are not all the same size.

Is there a way to automate this process at all, or is doing it manually my only option?

Fishnet option doesn't work, and I really can't find anything else to try. Each polygon would be split at the midpoint so I'm thinking there must be some way to save a little time.

(edited for clearer picture)


texts: How to iterate throught fields names in a table and change row values?

texts:
I got table like that one (Example).

I need to change 0 to 1 and -1 to 0 in each field except first and second.

The problem is that every time i have different number of fields and field name is changed too, but i can use wildcard    "FID_"   to select a field name.

I read some stuff about, should be like that:

 def test(field.name):
  fieldList = ListFields(%path%)
  for field in fieldList:
    if (field.name)==0:
      EQ=1
    else:
      EQ=0
    return EQ

%path% should be the value taken from model-builder (if that is possilbe).


texts:
I want to make a buffer from the river line with only the points that have over 20 m (elevation) from river line. (I hope you understand:D).
How can I do that?

Here is the picture....I want only the consecutive points over the river line :)

I use ArcMap 9.3. I've made offset points (from 50m to 50m) of river line from 200m to 200m. I interpolate with DEM and I query the points with 20m over the river line elevation...But I have a problem: I have a lot of points with 20m over, and I want only the consecutives one. How can I do this?

All the points with 20m over the river line
---
The points I want


texts:
I have some data in  PostgreSQL which I am exporting to Amazon RDS service.
Is it possible to use PostGIS extension in AWS?
(And will I need the mandatory geom column or I could run queries with only latitude and longitude values?)


texts: Can't upload .osc to local Rails Port instance

texts:
I have been experimenting with a local instance of Rails Port. I used osmosis to load a small pbf of my neighborhood. Using JOSM I edited some data on the local instance. Then a few weeks later I downloaded the latest osm of my neighborhood. Using osmosis I generated a .osc file, a diff between my postgres database and the new osm file. I opened the .osc file with JOSM where I could see the changes, but I couldn't upload the changes to local postgres database. It shows a conflict and upon giving sync, both entire data set and specific node it wouldn't upload ans keep showing the same conflict message.

Error message from JOSM: 'Conflicts detected Uploading failed because the server has a newer version of one of your nodes, ways, or relations. The conflict is caused by the node with id 1,472,108,803, the server has version 1, your version is 2. Click to Synchronize node 1,472,108,803 only to synchronize the conflicting primitive only. Click to Synchronize entire dataset to synchronize the entire local dataset with the server.Click Cancel to abort and continue editing.' If I click on any of them it syncs with the database and gives me the same conflict error.

I have a user account on local Rails Port instance with what I can upload data manually and I changed the OSM Server URL to point my local instance which is http://localhost:3000/api


texts: Select By Attribute: If-Then VB Statement

texts:
I am attempting to write an If-Then Statement where if all the UtilityBillingMerge101719.() are true then return all of the results that are true that have an account number.

This currently works, but am struggling to get it into an if-then format.


texts: How to select all remaining features after a Layer Definition filtering?

texts:
I have a map service that represents the World map with a table with population associated. I would like to filter some of the data matching a criteria and then, select ALL remaining features to work with them later, for instance, adding all them to an array.

For example, I can solve the first part of the problem if I create an ArcGIS Dynamic Map Service and apply a Layer Definition filtering to obtain a subset, like filtering countries with an population bigger/smaller than a threshold.

But, once I see on the screen my filtered map, how can I select all remaining features?

Does anyone know? Any Feature Layer magic?

Thanks for your help!


texts: File Browsing Dialog for Tables in Spatial Databases

texts:
I am looking for a file browsing object that allows a user to navigate into the different types of geodatabases.  For example, the ability to navigate into a Spatialite/PostGIS/ESRI File Geodatabase...etc database to select a specific spatial table.  I am working on a .NET based app.

In ArcObjects, the equivalent object would be the IGxDialog.


texts: I was wondering if any of you out there have experience with any apps that work on windows 8 tablets free or paid?We have a client that wants us to use their tablets for data collection.  i was wondering if some of you out there have used apps for the windows store that you can recommend.


texts: How can I sucessfully import 3D surface and features from ArcGIS 10 into Sketchup 8?

texts:
What are the steps required to successfully export an entire 3D model (composed of TIN surfaces, 3D polygons, 3D polylines and 3D points) from ArcGIS 10 to Sketchup 8?   Preferably, I would like the various components to line up. (Does not need to retain spatial reference)

Gone are the days when I was able to export TIN to Sketchup and easily export polygons, polylines and points with the ArcGIS plugin.

ESRI seems to think the process is a lot easier now... "Up until now this was a pretty cumbersome process. At ArcGIS10, this work-flow becomes much easier."

I respectfully disagree; Up until ArcGIS 9.3.1 it was a simple point and click operation (using the plugin) to select surfaces and whatever layers or selected features in the plugin and export. In fact it even opened SketchUp and loaded everything in and smoothed the surface if the check-mark was ticked off. And spatial reference was retained!

Has anyone tried to export complex TIN surface to Sketchup with the Collada technique? (No direct TIN export functionality now)
The Interpolate Polygon to Multipatch using TIN as input surface yields a result that looks nothing like the input surface:

As for line and point features it is possible to export to DXF 2007 then import into sketchup. Of course the data is not referenced so every imported DXF is placed in the default 0,0,0 position. I could overcome this by exporting "bounding box" features with each DXF but still - cumbersome....

I know people (myself included) have been asking similar questions but I still don't seem to be able to do this. Does anyone have a solution?


texts: OpenWind Zone of Visual Impact transparent pixels

texts: How do you make pixels which are not visible from a point be made transparent and the pixels which are visible opaque so that I can overlay the raster over a map to see clearly where you can and cannot see a wind farm from?


texts:
I am having a problem getting the results back from a WMSGetFeatureInfo call to my GeoServer. Using OpenLayers 2.12, GeoServer 2.4.0, Windows Server2008, and Microsoft IIS 7.5 with Application Request Routing Cache and URL Rewrite. I do NOT use a proxy because I use URL Rewrite. The URL Rewrite takes my incoming requests and formats them to avoid the cross domain/same origin problem.

I have been successfully querying WFS data using this setup for 2 years but when I try to use the WMSGetFeatureInfo it fails.

Here is my process:

I serve multiple WMS layers from GeoServer to my web page.
The user selects a WMS layer from a dropdown. This sets the variable "selLayer" so my WMSGetFeatureInfo call is pointing at the layer the user wants to get information about.
When the user clicks on the layer it fires this code which creates a OpenLayers WMSGetFeatureInfo control with an event listener that should accept the return text & xmlHttp request and display them in a div.  :

     var infoControl = new OpenLayers.Control.WMSGetFeatureInfo({
         url: 'http://xx.xx.xx.xx:8080/geoserver/wms',
         title: 'Identify features by clicking',
         layers: [selLayer],
         queryVisible: true,
         //eventListeners: { getfeatureinfo: function (event) { alert("TEXT: " + event.text + ", XY: " + event.xy + ", REQUEST" + event.request + ", FEATURES: " + event.features); } }
         eventListeners: { getfeatureinfo: function (event) {
             document.getElementById("divAttributes").innerHTML = event.text + "Test Text";
             alert("PUT TEXT HERE!!!");
         }
         }
     });

     map.addControl(infoControl);
     infoControl.activate();

But, when the script runs I do not get return values for the event.text. When I run this in FireFox (with FireBug) I can copy the GetFeatureInfo request into a new browser window and I get the desired return values.

So my question is: Why am I not getting the return values when I run this on my map webpage but I do get results by posting the request into a new browser window.


texts: How can I execute my custom model's tools in VB.Net?

texts:
I want to execute one of my tools (which I created with Model Builder).
I can execute a built-in tool by this code, but when I change the tool name with Model Tool Names it doesn't work.

I gathered information about referencing a new *.tbx file into a VB.Net project but it doesn't work properly. (When I use 'ArcGIS Toolbar Reference', nothing happens.)
Please help me.

This is my code for built-in tools:

Public Function CreateTempFileGDB() As Boolean

    Dim GP As ESRI.ArcGIS.Geoprocessor.Geoprocessor = Nothing
    GP = New ESRI.ArcGIS.Geoprocessor.Geoprocessor()

    GP.OverwriteOutput = True

    Dim parameters As IVariantArray = New VarArray
    parameters.Add(tempFolderPath)
    parameters.Add("tempFGDB")

    Try
        GP.Execute("CreateFileGDB_management", parameters, Nothing)
    Catch ex As Exception
        Return False
    End Try

    Return True

End Function


texts: Test the relationship between two point patterns

texts:
I have two point patterns obtained from medical images. Points are used to represent cell markers.
By visualization, the point pattern A is extremely possible to be a subset of point pattern B. So is there any statistical package or function in R that I can use to verify my guess?


texts: Creating seperate layers out of one big polygon shape by using attributes

texts:
I try to explain may problem. I have a big polygon shape (> 200 fields). I have the owners of these fields in a column. I need the seperated shapes for each owner. Is there a tool to do that?

Is it also possible to get an information about the whole area? Like a summary:
201 polygones
area: 1450 ha

In Saga I have a tool "Split shapes by attribute". This works fine with the expected result (10 owners: 10 new layers). This tool doesnt work with the sextante plugin.

Thanx!


texts:
I have a PostgreSQL database with a one-to-many relationship. Here is a reduced example:

CREATE TABLE buildings.building
(
  objectid serial NOT NULL,
  building_no character varying(40) NOT NULL,
  geometry geometry(Polygon,32638),
  CONSTRAINT building_pkey PRIMARY KEY (objectid)
)

CREATE TABLE buildings.flat
(
  objectid serial NOT NULL,
  fk_property_group integer,
  fk_building integer NOT NULL,
  CONSTRAINT flat_pkey PRIMARY KEY (objectid),
  CONSTRAINT flat_fk_building_fkey FOREIGN KEY (fk_building)
  REFERENCES buildings.building (objectid) MATCH SIMPLE
 )

I defined a relationship in the QGIS (2.2 and 2.4) Project Properties and it looks good in the forms. But when I enter a new building I will have the value
nextval('real_estate.property_group_objectid_seq'::regclass) for the objectid of the building
and the flats that I want to add won't have the correct foreign key. Its actually NULL.

Is there any workaround/solution for this issue?


texts: Measurement tool of QGIS working in unsaved project but disabled in saved project

texts:
The measurement tool of QGIS 3.8 is working in the unsaved project but disabled in the saved project. I unable change the measurements or the units, and the Units dialog box is grayed out.
However, when we open a new project the measurement tool is working as expected.
Why is this?


texts:
I would like to perform routing in MS SQL 2012 similar to pgRouting.

In pgRouting, I'd do the following steps to create a routable network:

Create a database routing with template template_routing
Create table road_network with following constraints

CHECK (st_ndims(the_geom) = 2)
CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL)
CHECK (st_srid(the_geom) = 4030)

Then
-- @ CREATE INDEX FOR THE ROAD TABLE -------------------- IMPORTANT

CREATE INDEX spatialindex_road
    ON road_network
    USING gist
    (the_geom);

Perform the following queries:

ALTER TABLE road_network ADD COLUMN "source" integer;
ALTER TABLE road_network ADD COLUMN "target" integer;

SELECT assign_vertex_id('road_network', 0.00001, 'the_geom', 'gid');

CREATE INDEX source_idx ON road_network("source");
CREATE INDEX target_idx ON road_network("target");

ALTER TABLE road_network  ADD COLUMN length double precision;
UPDATE road_network  SET length = length(the_geom);

ALTER TABLE road_network  ADD COLUMN reverse_cost double precision;
UPDATE road_network  SET reverse_cost = length;

ALTER TABLE road_network  ADD COLUMN x1 double precision;
ALTER TABLE road_network  ADD COLUMN y1 double precision;
ALTER TABLE road_network  ADD COLUMN x2 double precision;
ALTER TABLE road_network  ADD COLUMN y2 double precision;

UPDATE road_network  SET x1 = x(ST_PointN(the_geom, 1));
UPDATE road_network  SET y1 = y(ST_PointN(the_geom, 1));
UPDATE road_network  SET x2 = x(ST_PointN(the_geom, ST_NumPoints(the_geom)));
UPDATE road_network  SET y2 = y(ST_PointN(the_geom, ST_NumPoints(the_geom)));
alter table road_network add column cost double precision default 0;
update road_network set cost=0.1 where type='NH';
update road_network set cost=0.2 where type='SH';
update road_network set cost=0.3 where type='major';
update road_network set cost=0.4 where type='minor';
update road_network set cost=1.2 where type='colony';
update road_network set cost=0.8 where type='third';

Now the network table is created.

To check, run this query on your postgreSQL....
assign table name name and create a table

CREATE TABLE shortest_path_astar_table_3(gid int4) with oids;
SELECT AddGeometryColumn( 'shortest_path_astar_table_3', 'the_geom', 4030, 'MULTILINESTRING', 2 );
INSERT INTO shortest_path_astar_table_3(the_geom)
SELECT the_geom FROM astar_sp_directed('road_network',37,43,true,true);

Open shortest_path_astar_table_3 on QGIS and check the path.

Is there any similar way to perform similar queries in SQL Server 2012?


texts: How to detect whether any points in the coastline is facing to the open sea or not?

texts:
I have a set of Mangrove stands polygon derived from remotely sensed data and the coastline data (Line). I am going to detect any point in the boundary of mangrove polygon is located in the outer side of the coastline and exposed to the open sea?

Can anyone suggested me how can i do it in ARCGIS?


texts:
I am trying to run:

SpatialReference poSourceSRS = new SpatialReference();
poSourceSRS.ImportFromEPSG(4326);
String unit = poSourceSRS.GetAttrValue("UNIT");
System.out.println("unit: " + unit);

on a Tomcat 8(Run As.. -> Run on Server), which I linked to eclipse luna, and I get the tomcat error(in eclipse):

SEVERE: Servlet.service() for servlet [package.Servlet] in context with path [/myServlet] threw exception [Servlet execution threw an exception] with root cause
java.lang.ClassNotFoundException: org.gdal.osr.SpatialReference
    at org.apache.catalina.loader.WebappClassLoaderBase.loadClass(WebappClassLoaderBase.java:1305)
    at org.apache.catalina.loader.WebappClassLoaderBase.loadClass(WebappClassLoaderBase.java:1157)
    at package.Servlet.showTileURL(Servlet.java:45)

But when I try to do the same in a java application (Run As... -> Java Application) it works..

I read something about editing catalina.bat but they used Linux and I am on a Windows 7 64-bit. I don't know what I have to add to catalina.bat


texts: MapInfo to layered pdf with linked layered legends

texts: I am trying to make a layered pdf of Rock geochemistry thematic map from MapInfo of several different elements on same base so I can present a clean 1 page pdf map with the ability to view or unview the different elements values with the appropriate legend with one click.  I get it all to work except the legends will not come in as separate layers attached to proper map. Can anyone help?


texts: Changing colours in SVG symbol using QGIS?

texts:
I am using QGIS version 2.14.1 (Essen), and I have imported a traffic cone symbol in SVG format to represent some of my data. I was wondering if it was possible to change the color of this symbol according to the values in one of the columns of my shapefile?

Basically my traffic cone is orange and white, and I would like to change the shade of the orange according to the values in a column of my shapefile. Is it possible to do this in QGIS?


texts: How to use filter in GeoExt.FeatureStore?

texts:
I have a FeatureStore:

        var protocol=new OpenLayers.Protocol.HTTP({
            url: "/geoserver/ows?service=WFS&request=GetFeature&typeName=xxxx:lineview&srsName=EPSG:4326&outputFormat=GML2",
    });
    zoomStore_disp = new GeoExt.data.FeatureStore({
        layer: myVecLayer_disp,
        fields: [
            {name: 'id', type: 'int'},
            {name: 'name_dispatcher', type: 'String'},
            {name: 'voltage', type: 'String'},
        ],
        proxy: new GeoExt.data.ProtocolProxy({
            protocol: protocol,
        }),
        autoLoad: true,
    });

But lineview layer contains too much features. So I must wait a lot of time before store is load. Is it possible to load in store features by filter?

I can't find any examples or documentations how to use it.


texts: OpenLayers 2, GeoJSON can't get coordinates to work

texts:
I'm trying to display some polygons on a map, and I'm new to OpenLayers, maps and GeoJSON.

First what works:

vector1 = new OpenLayers.Layer.Vector("GeoJSON1",
{
    projection       : "EPSG:4326",
    strategies       : [new OpenLayers.Strategy.Fixed()],
    protocol         :  new OpenLayers.Protocol.HTTP({
                            url: "test.php",
                            format: new OpenLayers.Format.GeoJSON()
                        })
});

This piece of code works, the GeoJSON from test.php shows up on the server, at the position where it should be. But GeoJSON I try to reach is on another server and it throw errors (I know I can set headers and it would work) but I don't want to do it this way.

This doesn't work:

var geojs_format = new OpenLayers.Format.GeoJSON();
var geojsval= {
                "type": "Polygon",
                "coordinates":
                    [[[11.0878902207, 45.1602390564],
                      [14.931640625, 40.9228515625],
                      [0.8251953125, 41.0986328125],
                      [7.63671875, 48.96484375],
                      [11.0878902207, 45.1602390564]]]
            };
vector = new OpenLayers.Layer.Vector("GeoJSON",{projection:"EPSG:4326"});
vector.addFeatures(geojs_format.read(geojsval));

This does not work - It displays polygon at a different place even though projection is specified in same way (or at least for me...).


texts: QGIS snaps with offset while using advanced digitization

texts:
I've ran into a problem, where if I draw two parallel lines, enable Advanced Digitization and start a perpendicular line from one line to the other with Snapping enabled, when I try snap to the intersection, the perpendicular line moves up or down with an offset.

It snaps fine if digitization is not enabled, or if I want to snap on any other location, however as soon as my mouse gets to the expected intersection point, the target shifts.

Edit 1:
I found another person having a similar issue here


texts:
I have a .GPB file, I need to extract the common GPS data from it but I don't find a clear structure that specify the contents of that file.

as a sample record from file

00000260 02 11 A9 61 E7 96 C4 95 B6 86 CA 01 9A 09 DC 01

how can i go through this file to get my data, I searched alot and found that it called "Dell Axim Navigation System (.gpb) file format".

Regards,

m.sai


texts: Publish a Locator Service with a place-name alias table

texts:
I am trying to publish a locator service with support for a place-name alias table, in ArcGIS for Server 10.3.

It should function similar to the Esri World Geocode Server which allows users to search for place names such as "Starbucks" over an extent.

So far I have not been able to publish the locator when the alias table is attached.  I found this article which describes how the place-name tables are not supported with published services.

Does anyone know how Esri was able to accomplish this with the World Geocode Server?


texts: Paint a FeatureCollection using color values stored in a property of the FeatureCollection itself

texts:
Starting with an example shown here: https://groups.google.com/g/google-earth-engine-developers/c/w20Unaq6zCw/m/a8HrcpqSCAAJ
How would I modify the code below to Paint using a Color value stored as a property in the FeatureCollection.
So, for example, in the regions variable below, {name: 'Rice'} is instead stored as {Color: 'blue'}
Then, rather than having the palette variable hard-coded as below, the palette would come from this Color property stored in the feature.  Note that I won't know ahead of time how many different colors are stored in the featurecollection.
This seems like it should be straightforward, but I cannot get it to work.
//specify regions corresponding to different vegetation types, here as representative rectangles with specified coordinates
var regions = ee.FeatureCollection([
 ee.Feature(ee.Geometry.Rectangle(-121.657, 38.109, -121.654, 38.108), {name: 'Rice'}),
 ee.Feature(ee.Geometry.Rectangle(-122.024, 38.199, -122.02, 38.197), {name: 'Annual grassland'}),
 ee.Feature(ee.Geometry.Rectangle(-121.875, 37.847, -121.873, 37.8455), {name: 'Forest'}),
 ee.Feature(ee.Geometry.Rectangle(-121.803, 37.8615, -121.8015, 37.8603), {name: 'Shrubland'}),
 ee.Feature(ee.Geometry.Rectangle(-121.825, 38.049, -121.823, 38.0473), {name: 'Wetland'}),
 ee.Feature(ee.Geometry.Rectangle(-121.825, 38.049, -121.823, 38.0473), {name: 'Wetland'}),
 ee.Feature(ee.Geometry.Rectangle(-121.825, 38.049, -121.823, 38.0473), {name: 'Wetland'}),
]);

function categorize(fc, prop) {
  var uniqueValues = ee.List(fc.distinct(ee.SelectorSet(prop)).sort(prop).aggregate_array(prop))

  return fc.map(function(f) {
    return f.set(prop + '_category', uniqueValues.indexOf(f.get(prop)))
  })
}

var count = regions.aggregate_count_distinct('name').getInfo();

var palette=['FF0000','00FF00','0000FF', '8b4500', '000000', '00FFFF'];
var regionsimage=ee.Image().toByte();
var paintimage=regionsimage.paint(categorize(regions, 'name'), 'name_category', 4);
Map.addLayer(paintimage, {palette:palette, min:0, max:count});

Map.centerObject(regions)


texts: Detecting click on map outside features in OpenLayers 2?

texts:
How to detect click on map outside features from two vector layers in OpenLayers 2?
Here is a code that alerts on selecting feature:

var layer1 = new OpenLayers.Layer.Vector("Vector Layer 1");
var layer2 = new OpenLayers.Layer.Vector("Vector Layer 2", {
    styleMap: new OpenLayers.StyleMap({ "default": {
    'fillColor': '#669933',
    'pointRadius': 6
  }})
});
for (var i=0; i<10; i++) {
    layer1.addFeatures([new OpenLayers.Feature.Vector(
    new OpenLayers.Geometry.Point((360 * Math.random()) - 180, (180 * Math.random()) - 90)
  )]);
    layer2.addFeatures([new OpenLayers.Feature.Vector(
    new OpenLayers.Geometry.Point((360 * Math.random()) - 180, (180 * Math.random()) - 90)
  )]);
}
map.addLayers([layer1, layer2]);

var select = new OpenLayers.Control.SelectFeature([layer1, layer2], {
    onSelect: function (feature) {
    alert("selected");
  }
});
map.addControl(select);
select.activate();

Now I want to alert "all deselected" when user resets the selection.

Example:
http://jsfiddle.net/f6rvj7q9/


texts: Identifying cluster centroids using ArcGIS Desktop?

texts:
I have a data set of points which are grouped into clusters like so:

I'd like to use ArcGIS or Python to create a point layer of the centroids of each cluster but I can't seem to find the right tools. I looked into hot spot analysis but it doesn't seem to apply, since my features aren't numerically weighted. I also tried the Integrate tool with a reasonable tolerance but I didn't see any changes.


texts: Editing existing polygon vertices to create smooth curves using ArcMap

texts:
I'm editing a really detailed polygon in ArcMap. The vertexes of the polygon are connected with an straight line. I want to convert them in a smooth curve.
While creating the polygon I use the "End Point Arc Segment" tool to create smooth curves:

But for the already created polygons, this tool is disabled.
How do I use this tool to correct the imperfection in my polygon?


texts: Working with Post miles and aerial imagery

texts:
As a highway design engineer most of the time we work with county routes and post miles.

For persons who don’t know what a post mile is it is basically a mile marker number on the highway of where you are. These numbers increase from south to north and west to east.

For example, I would have something like this Mendocino 20 PM 7.5 – 14.5 with this, I would like to know if qgis can search for me an aerial imagery within the area.

I am not good in writing programs so something more automatic maybe with a plugin that has been made already would be great.


texts: Format Query statement for use in Virtual Layer query from PyQGIS

texts:
Based on QGIS: SQL query within python console to sum a column based on values of other two columns

I am trying to write a script that will get the field names to use from PyQGIS but am having issues with the format of the query statement.

I have

layer = qgis.utils.iface.activeLayer()
for field in layer.fields():
    current_field=field.name()
    print(current_field)
    query = "Select", current_field, "count(*) from", layer,"group by", current_field,"order by 2 Desc"""
    print(query)
    vlayer = QgsVectorLayer( "?query={}".format(query), current_field, "virtual" )
    QgsProject.instance().addMapLayer(vlayer)

But this is not formatting the query correctly so it comes out as

('Select', 'Type', 'count(*) from', <qgis._core.QgsVectorLayer object at 0x00000244813FD558>, 'group by', 'Type', 'order by 2 Desc')

Which gives the following error

Query execution error on CREATE TEMP VIEW _tview AS ('Select', 'Type', 'count(*) from', <qgis._core.QgsVectorLayer object at 0x00000244813FD558>, 'group by', 'Type', 'order by 2 Desc'): 1 - near "(": syntax error

How do I format the query= correctly? I have also tried it with layer.name()

It should be formatted as

Select Type, count(*) from SW_Test group by Type order by 2 Desc

Which results in


texts: Export Feature Class to Geodatabase without SubTypes

texts:
I want to export a feature class to a geodatabase using arcpy and strip out any subtypes.  But they keep hanging around in the new feature class.  How do I get rid of them?  This doesn't work:

subtypes = arcpy.da.ListSubtypes(featureClass)
arcpy.RemoveSubtype_management(featureClass, subtypes)

The above throws this error:

Traceback (most recent call last):
File "E:\SCAT\Scripts\test.py", line 112, in <module>
arcpy.RemoveSubtype_management(featureClass, subtypes)
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 17276, in RemoveSubtype
raise e
RuntimeError: Object: Error in executing tool


texts:
I am building a composition with map, legend... I want to add a paragraph (text) to the composition but the text only recognize one word. What I am doing wrong?

Code:

text = QgsComposerLabel(myComposition)
text.setText("Prognostico" + "\n" + "Coordinate System: ETRS89 PTTM06" + '\n' + str(self.pr))
text.setFont(QFont("Cambria",15, QFont.Bold))
text.setItemPosition(328.2, 72.3, 87.7,102.3)
text.adjustSizeToText()
#title.setFrameEnabled(True)
myComposition.addItem(text)

In this example only Prognostico word is showed in the layout. Why the rest doesn't appears?


texts: Can I see what hillshade illumination settings were used to create a particular hillshade geotiff?

texts:
Are the settings used to create a hillshade geotiff captured in the geotiff file somewhere? I am trying to expand the coverage of an existing hillshade and want to make sure I am using the same Illumination settings (Azimuth, Vertical Angle) as the existing portion.
Any idea where I could look up that info?


texts:
I have a layer containing linestrings and want to move those that are too close to each other.

I seem to be missing something in the moving part though:

layer = iface.activeLayer()

def moveFeature (fId,dx,dy):
    layer.startEditing()
    layerUtil = QgsVectorLayerEditUtils( layer )
    result = layerUtil.translateFeature(fId,dx,dy)
    layer.commitChanges()
    layer.triggerRepaint()
    return result

since the coordinates of the first point of any feature before and after move are the same (thoug the feature is moved in the map view) ie the piece of code:

for feature in layer.getFeatures():
    if feature.attribute("IDR_LIT") == "EGUZOL71RUEYR":
        print feature.attribute("IDR_LIT")
        print "before move coordinates: %d;%d"%(feature.geometry().asPolyline()[0][0],feature.geometry().asPolyline()[0][1])
        moveFeature(feature.id(), 1111,1111)
        layer.updateFeature(feature)
        print "after coordinates: %d;%d"%(feature.geometry().asPolyline()[0][0],feature.geometry().asPolyline()[0][1])

prints:

EGUZOL71RUEYR
before move coordinates: 675392;6408936
after coordinates: 675392;6408936

Does anyone have a hint about this?


texts: Unable to preventDefault inside passive event listener due to target being treated as passive: Openlayers 2

texts:
I am getting the following error messages: How do I define add control and prevent passive event listener?

Below is the code:

<!DOCTYPE html>

 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>

<button dojoType="dijit.form.ToggleButton" id="featureInfoButton"
onChange="featureInfoChange" iconClass='dijitCheckBoxIcon'
checked="false">Activated</button>

Additional question is regarding this line:

 function init() {
 map = new OpenLayers.Map('ch05_wmsfeatureinfo', {});

what is ch05_wmsfeatureinfo, is this the wms layer one trying to query? I am following the book tutorials. Please see book code at this link: No GetFeatureInfo response when clicking on a map in OpenLayers 2


texts: Which EPSG for UTM with 8 digits in QGIS

texts:
I have a UTM datasets with 8 digits, so the first digits are 32 (for UTM Zone 32N), but which EPSG should I use to get the data at the right spot in QGIS in a 6 digits system.

If I want to transfer UTM coordinates with 6 digits to 8, I think I could use EPSG 5652 (ETRS89 / UTM zone 32N (N-zE)) but which one I should use to get a 8 digit value to the 6 digits type?


texts:
I have two vector line datasets and wish to compare the line lengths. While I have already compared the sum of the lengths, I'd like to conduct further analysis down to 250m, possibly lower.

I'm using qGIS, and already know how to create a vector grid. However I would like to create a grid which clips per individual square and compares the length of the two datasets contained within.

I realise that this is a process that could be completed manually, however the size and complexity of the data (all roads at a city level) mean manually doing the task would take months.


texts: Creating symbols between connecting lines to represent line boundaries

texts:
I am trying to figure out a way to create markers as seen from the attached image where an angled line (in dark red) is created between the lines to represent a connecting stream. With the data I am using, each break is actually an assessment unit and all assessment units form a segment (stream). Does anyone know how these are created? I would assume an "end of line" symbol" but not sure where to look for that.

I am using ArcMap 10.6 and have access to most all tools.


texts: Finding points along line closest to known point

texts:
I'm working on a TV graphics system which will be used on a broadcast from train ride.
With this system I want to automatically be able to list interesting POIs along our route.
We have a large postgis table of POIs which are within 2km of the train route, and we have a GPS fitted to the train, so we know both our current position, the polyline of the train route and the lat/lng of the POIs close to the route.
I've been wrestling with understanding how I can most easily get a hold of the POIs which are closest to my current posistion at any given time.
I have two PostGIS tables:
The first is called 'tracks' with two columns:
date (date) | track (polyline)

The other is called 'pois', and has several columns, but the relevant one is
the_geom

I have a GPS connected by socket which I can read whenever.
So the query I'm interested in formulating is something like:
"Give me the POIs closest to my current GPS position along this linestring"
I've been looking at using ST_BUFFER in combination with ST_DWITHIN to get it up and running, but I'm obviously not at a high enough competency level to get it working. Seeing as both my GPS position and the POIs position probably will never be exactly on top of the linestring, will this maybe complicate things further?
Has anyone done something like this?


texts:
I can't seem to figure out a way to define a custom PROJ4 projection that can mirror geometries (e.g. inverse x to -x and y to -y and vice versa). Other than that, I don't need any transformations or reprojections. I am simple trying to show arbitrary (x,y) coordinate pairs from an SVG in a GIS system (OpenLayers and/or QGIS).

I think the main problem is that in SVG, the origin (0,0) is situated at the Top-Left. However, I have no clue how to account for that difference..

I have seen that PROJ4 allows transformations (+towgs84) but I figured that's not what I should use. I know from visualizing my geometry in QGIS and OpenLayers that a 'mirror functionality' is all I need.

I am open for other solutions than projections, but I figured that defining a new projection would be easiest. I do need an automated way though.

TL;DR:
I would like to mirror-transform (projection or something else) X,Y-coordinates for using in QGIS or OpenLayers without use of external libraries or plugins.


texts: We need to flip a some vector geometry, like a pancake. We've tried in QGIS 3.10 the Shape Tool, but it's not working. Do any of you have an idea or know a plug-in that would work?


texts: Seeking LAG LiDAR Analysis GUI installation guide?

texts:
Does anyone know how do I install LAG lidar Analysis Gui?
LAG by arsf

This is my situation: I'm using debian 8 test. I installed this library: boost (yes) * GTKmm (yes) * GTKGLextmm (yes) * GThread (yes) * libgeotiff (no) * laslib (yes with Lastool) * lidarquadtree (no). After I Tried this: 1) Install the binutils-gold package before compiling but it returned an error (dependency not...). I stopped here.


texts: Integrating application security with GeoServer security?

texts:
I'm going to develop a map server with my own logic and entities. It will have a Postgres database, user management, specific layers with certain types, WFS, WMS, etc.
I'm going to use springframework and GeoServer as the GeoServer is an open source project.
I have user a management problem. How can I integrate my own project security with that of GeoServer? I mean how can I add a user or role or group or assgin a role to a user programmatically? How can I authenticate with geoserver programmatically?
As you know you can manage layers, workspaces, etc by GeoServer-manager. But it does not have user management REST API!

I've developed a map server that maintains layers info and I want it to support WFS and WMS. I want to grant access to my users of layers.
It was the most details of my application.


texts: Plot points on map using Fiona and/or Shapely

texts:
I have written a parser in Python that extracts the coordinates out of a GPX file and I would like to view those points on a map.  I have done a little research and it seems like either Fiona or Shapely is capable of doing this, but all of the examples I find handle shape files.

My coordinates are stored as tuples in a list.  I was wondering if there was an easy way to use either of these packages to plot my coordinates.


texts: Interpreting variogram: negative values beyond model sill?

texts:
The negative values mean that there is negative autocorrelation at some far distance, like a repulsion.  How would this affect assessment of the range?

(This is just for learning the scale of a spatial process, not for kriging.)


texts: QGIS Offline editing of PostGIS table produces a Spatialite table with no attribute lengths?

texts:
When synchronizing an offline Spatialite file with it's corresponding PostGIS database, using the Offline Editing plugin in QGIS, I received an error I hadn't experienced before

The error message is quite clear, I had entered more than 254 characters in a attribute field that in the PostGIS table had this field length.

I looked at the properties of the offline Spatialite table and all the attribute length and precisions are all '0', or undefined.

In the original PostGIS database table they are defined

Even if the Spatialite attribute that is too long is corrected, the synchronising of the offline Spatialite file now gives a commit error, which means all the changes made offline will now not synchronise.

My question

Is it possible to bring across the attribute field lengths from the PostGIS table to the Spatialite file so as to not allow attributes that are too long to be entered in the first instance?


texts: Formula for Nicolosi Globular Projection

texts:
What is the formula for the Nicolosi Globular projection?

In other words, what is the function that maps the coordinate (φ, λ) on the sphere to the coordinate (x, y) on the plane?

Example of a Nicolosi Globular projection

You can assume that there are two such functions, one for each of the two hemispheres of the projection, such that the origin (x = 0, y = 0) in both cases is in the center of the given hemisphere.


texts:
I read a .csv file with geometry like this:
trajs = pd.read_csv("trajs.csv")
trajs["geom"] = trajs["geom"].apply(wkt.loads)
trajs = gpd.GeoDataFrame(trajs, geometry="geom")
ax = trajs.plot(figsize=(10, 10))

But it shows like this:

But in QGIS, it show like this:

How can I set the range of GeoPands plot?


texts:
I am trying create a topology for a zoning dataset in my SDE feature dataset. I've run into a few issues. First off was an issue with versioning that this questions resolved: How to solve "Cannot create topology" in ArcCatalog?

Now I've run into a second (well more like fifth) issue concerning the spatial references of the Topology / Zoning dataset:

The feature dataset has the pre-defined projection we use for all our data. In terms of spatial references, all layers within the feature dataset are governed by the spatial reference / projection of the feature dataset correct?


texts: PostGIS where clause not working as expected with multipart geometry

texts:
I am running an intersect from a subset of one table on itself (this is part of a larger parameterized query but the issue in question arose out of this particular situation) using a CTE (zones in the query below).
zones should result in four records selected by ID.
In the subsequent query summ_table should be limited to just three of those four records due to the where clause (WHERE (ST_Area(summ_table.geom) > 7769967522)).
However, when I run the query it produces the following result:
id  area
4   7911449533.18243
36  8167237869.65092
38  5.16730218436123e-05
69  17451037319.6371

id=38 is obviously not correct. I assumed the result would only include the three records that satisfy the above where clause.
Here is the full query:
-- select from zone table using IDs
WITH zones AS
(
    SELECT ogc_fid, geom
    FROM mn_county_dissolved
    WHERE ogc_fid IN (4,36,69,38)
)

-- individual area of features that intersect each zone
SELECT inter.ogc_fid, sum(inter.area) AS area
FROM
(
    -- zone / summary table intersect
    SELECT
        zones.ogc_fid,
        CASE
            WHEN ST_CoveredBy(summ_table.geom, zones.geom) THEN ST_Area(summ_table.geom)
            ELSE ST_Area( ST_Intersection(summ_table.geom, zones.geom) )
        END AS area
    FROM mn_county_dissolved AS summ_table
    INNER JOIN zones
    ON ( ST_Intersects(summ_table.geom, zones.geom) )
    WHERE (ST_Area(summ_table.geom) > 7769967522)
) AS inter
GROUP BY inter.ogc_fid

And here are the values for those records in question for ST_Area:
 ogc_fid | countyname  |    ST_Area
---------+-------------+------------------
       4 | Beltrami    | 7911449533.18243
      36 | Koochiching | 8167237869.65092
      38 | Lake        | 5924950382.16346
      69 | St. Louis   | 17451037319.6371

I'll help you out with the math --> 7769967522 - 5924950382.16346 = 1845017139.84
As you can see the record with id=38 should definitely not satisfy that where clause
The only thing I can think of is that it might have something to do with the number of geometries for each record.
ST_NumGeometries returns 8 for the record id=38 and 1 for all of the others.
Note that ST_GeometryType returns ST_MultiPolygon for all rows.
Why is that where clause in the query not filtering out that particular row (id=38), and if the query is indeed working as expected, is it due to the fact that there is more than one geometry for that particular record?


texts: Python scripting for making rings from the city center

texts:
Problem description:

I have two layers: 1) [Polygon] city blocks, and 2) [Point] the city center. The City blocks layer has 110,681 records. And The city center layer has just one record because there is only one city center in this city.

arcpy.SelectLayerByLocation_management("CityBlocks","WITHIN_A_DISTANCE","CityCenter","1 Kilometers","ADD_TO_SELECTION")

The code above selects census blocks that are 1 kilometer away from the city center in ArcMap.

Question:

Now, I have to add a field called "ring" onto the city blocks layer. And, I have to write a code that does:

ring = 1 if city blocks are 1km away from the city center.
ring = 2 if city blocks are in the ring of 1km to 2km from the city center.
ring = 3 if city blocks are in the ring of 2km to 3km from the city center.
...
It stops if there are no city blocks in the ring of (N-1)km to (N)km from the center.

How do we write a code for that selection process using arcpy?

Thanks a ton in advance!


texts: Convert OSM File to get addresses using Osm2pgsql

texts:
I am working to covert OSM File to get all building's informations in the city.

I am using osm2pgsql package to transform the data. Now I am confused as I have 5 tables :

   - planet_osm_point
   - planet_osm_line
   - planet_osm_roads
   - planet_osm_polygon

https://github.com/openstreetmap/osm2pgsql

I am confused now, how can I get buildings' informations from these tables. What do every table represent ?


texts: What are the different algorithms employed for the generation of a Digital elevation model from images captured using an unmanned aerial vehicle ?


texts: How to read a heat map at a point in QGIS?

texts:
I have created a set of heatmaps based on the amount of certain chemicals released at specific sites in an area using the QGIS heatmap plugin. I now have a CSV of patient locations that I have overlayed onto this map.

I want to be able to read off the individual raster value (color/chemical concentration) for every patient at their exact location and be able to export this for analysis. How can I do this?

Additionally, I have heatmaps for almost 100 different chemicals so would it be possible to overlay all of the heatmaps and get an individual value for each chemical that could be appended to a new row on the patient attribute table? If this is too difficult or not possible that is fine.

Thank you so much for your help!


texts: When i load my KML on google earth, the balloon appears whole but when I embed on a webpage using the embedKML tool, the data in the balloon (In this case a table) appears to be cut. The close button is not showing. I have tried increasing width and height in the tool but there is no change.


texts: QGIS 3.4 using pygis for deleting colomns

texts:
I am currently building a model for importing and preparing multiple data by merging shapes and joining tables.

Almost finished with the model, but I haven't found a useful tool for deleting certain columns of an attribute table in the graphical modeler.

So I am right now stuck in writing myself a tool for this purpose and unfortunately I am not that experienced in writing code.

I just need to delete following indexes - Always the same ones.

eg. :[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,19,20,22,23,24,25,27,28,29,41,42,44,45]

edit: Code:

layer=Qgs.utils.iface.activeLayer()

layer_provider.deleteAttributes([2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,19,20,22,23,24,25,27,28,29,41,42,44,45])
layer.updateFields()

edit2: I tried to make a dropdown menu using an existing code from J. Monticolo It did not work, as expected..

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterFeatureSink
import processing

class MyModel2(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('inputlayer1', 'input_layer_1', types=[QgsProcessing.TypeVector], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Output_layer', 'output_layer', optional=True, type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
        results = {}
        outputs = {}

        layer=QgsVectorLayer('input_layer_1')
        with edit(layer)
            layer.deleteAttributes([2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,19,20,22,23,24,25,27,28,29,41,42,44,45])

        return results

    def name(self):
        return 'delete_table'

    def displayName(self):
        return 'delete_table'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return MyModel2()


texts: Create my own MapBox Vector tile and self hosting

texts:
I understood hosting mbtiles is possible and some solutions are available. I not sure that I can find the answer;
How can I create my own vector tiles (mapbox vector tile format) which from my data stored in postgis or geojson etc.

Should I download mbtiles from mapbox?


texts: Unable to Write to Text File Unless Python Script is Re Pointed

texts:
I have a python script running from ArcToolbox. Each time the script is run, it creates a new text file (testtext.txt), writes to it, then deletes it at the end.

Before the script is re-pointed, the source looks like this:

When I try to run the tool, I get the error:

If I go into the properties of the script and re-point the source to the exact same location as shown above, the script does work, and works fine every time thereafter. This problem happens when logged in as Administrator on the computer as well. Is there a way to ensure permission without having to re-point the source each time the tool is added to ArcMap?

I have "Store relative path names" enabled on the script properties. Here is the code snippet:

tempPath = sys.path[0]
newPath = os.path.join("testtext.txt")

open(newPath, 'w').close()
with open(newPath,"a") as newfile:
    newfile.write("http://genericaddress/#/selections/place?query=")
outfile = open(newPath,"a")


texts: Loading Landsat ImageCollection into GEE

texts:
Am trying to use this ee.ImageCollection() method to load Landsat data using Google Colab. See code below:

def maskL8sr(image):
  cloudShadowBitMask = 1<<3
  cloudBitMask = 1<<5
  qa = image.select('pixel_qa')
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).aAnd(qa.bitwiseAnd(cloudsBitMask).eq(0))
  return image.updateMask(mask)

l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2015-01-01', '2017-12-31').map(maskL8sr)
vizParams = {bands:['B4','B5','B2'], min: 0, max: 3000, gamma: 1.4,}
Map.setCenter(114.0079, -26.0765, 9)
Map.addLayer(l8sr.median(), vizParams)

But I keep getting this error:

Is it a bug?


texts:
I want to generate tiles for an PNG file. In the MapTiler(GUI wrapper of gdal2tile), I can specify the bounding box and the zoom level, then it does the good job, which generates the tiles only for that bounding box.

I want to do it with gdal2tiles with own control. How can I directly use gdal2tiles.py to specify the bounding box? Which other gdal utility programs should be called before gdal2tiles?


texts: Get invalid blocksize error when trying to generate a render base with osm2pgsql

texts:
I'm trying to generate a database with the osm2pgsql schema using the osm2pgsql application, with the command:

osm2pgsql -H host -U user -P 5433  -W -d db -k -l -c -s file.pbf

But when I'm running osm2pgsql, I get the following error:

Using built-in tag processing pipeline
Using projection SRS 4326 (Latlong)
Setting up table: planet_osm_point
Setting up table: planet_osm_line
Setting up table: planet_osm_polygon
Setting up table: planet_osm_roads
Allocating memory for dense node cache
Allocating dense node cache in one big chunk
Allocating memory for sparse node cache
Sharing dense sparse
Node-cache: cache=800MB, maxblocks=102400*8192, allocation method=11
Mid: pgsql, scale=10000000 cache=800
Setting up table: planet_osm_nodes
Setting up table: planet_osm_ways
Setting up table: planet_osm_rels
Reading in file: snapdb.pbf
Invalid blocksize 1010792557
Error occurred, cleaning up


texts: QGIS and ArcGIS common edit for Postgresql

texts:
I have created a database with postgis(v2.2.1) template in postgresql 9.5 version. Once created database, I used  Create Enterprise Geodatabase geoprocessing tool existing database via Arcgis (10.5).

And now I can create new feature in database with ArcCatalog. But I can not use this feature table in QGIS.

If I create a geometry table with plain sql like this:

create table mypoints (id serial, name varchar, geom geometry(Point, 4326));

this time arcigs will not use mypoints table.

I could not create a common solution. How can I do this? Or is tihs possible?


texts: How can I transform OSM road data to a graph network in 2018?

texts:
I see a lot of questions in this topic, but most answers are at least 4-5 years old. I wonder if there is a more current solution that would allow me to gain a graph network from the OSM highways data.

Simplicity is key.

Previous discussions:

How to Build a graph from a map

Read and Export OSM Road Data as a Network of Road Segments [Adjacency Matrices]


texts: FME Desktop saves Transformers uncommented in Workspace File

texts:
I'm using FME Desktop 2014 SP1 to create a rather huge mapping process from ArcSDE Feature classes to INSPIRE GML Writer. After saving the Workspace and reopening it a couple of days later, all Transformers and also some readers were gone.

I opened the workspace file with notepad++ where I could recognize some of the disappeared transformers information, but it was commented out with #

I also tried to save the workspace under different names, but every time the same happened.

I think it should be possible to recover the transformers by manually deleting '#'s but that's no suitable solution to keep on working on my workbench.

Had anybody issues like these or knows a solution to this?


texts: How to keep column headers with attribute table join

texts:
I successfully joined an attribute table to a shapefile using the method below.

But the new attribute table does not include the original column headers from the joined table. They've been replaced with something generic. See below.

Is there a way to keep the original column headers from the table that was joined?


texts:
I have 2 overlaping shapefiles. Is there any way, where I can select polygons from the 1st shapefile and trim everything that overlaps with them, in the 2nd shapefile?

I am using QGIS.


texts: I was going to add a field to my polygon layer for obtaining the area of the polygon and got an error message shown below. Never encountered this before - what is going on here?


texts: Concatenate Label Values Fails to Calculate Aggregate

texts:
I'm trying to do what seems like a very simple operation (It's SUPER easy in ArcMap...). I'm kind of a novice with QGIS, so it's probably something really simple that I'm doing incorrectly. The release is 2.18.12. Basically, I want to create concatenated labels with values from three separate fields, separated in the text by a comma and a space. The fields I'm trying to concatenate are numeric and have no null values (unless zero is considered null). Here is my expression:

Here is the error summary:

This is ridiculously easy in ArcMap using a JavaScript expression. Heeeeelp!


texts: Changing order of attribute data in InfoWindow on click of multiple features

texts:
I have an ESRI Web App Builer 2.7 stand alone web application where it is required that we are able to control the order of attribute data appearing in an Infowindow on clicking multiple features.
 For example, on click at a point on the map there might be several features overlapping each other. The ESRI out of the box infoWindow allows us to navigate through associated attribute data of each overlapping feature by click Next/Previous button on top right corner of the InfoWindow popup.

Is there any way to control the order in which the feature data appears? Please refer to attached image.


texts:
I have a point layer from an imported shapefile that I am trying to label automatically with arcpy. I had hoped that a simple:

def update_labels(id):
    import arcpy
    mxd = arcpy.mapping.MapDocument("C:\Users\username\Desktop\LOCATES.mxd")
    df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
    updateLayer = arcpy.mapping.ListLayers(mxd, id + '_fishnet_label', df)[0]
    sourceLayer = arcpy.mapping.Layer(r'C:\\Users\\username\\Desktop\\Before you dig\\label_template3.lyr')
    arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, False)
    arcpy.RefreshTOC()
    arcpy.RefreshActiveView()

would do something, but it doesn't seem to - no error but nothing happens. (template3.lyr is a point layer set up with the appropriate symbology and labels, all fields match. This example is run in IDLE with the arcpy module imported, but I have tried to run it in console with similar results)

Another suggestion was to change the data source of my template layer:

"label_template"

to a shapefile with point data:

"C:\\Users\\marcus.rodger\\Desktop\\Before you dig\\temp\\fishnet_label.shp"

using

arcpy.label_template.SetDataSource(mxd,r"C:\\Users\\username\\Desktop\\Before you dig\\temp\\fishnet_label.shp)

which just gives me a runtime error

I have been studying the help e.g. https://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-mapping/updatelayer.htm
but I'm not sure if I am just misunderstanding how it should be applied.

I am pretty new to arcpy and coding in general, the rest of the script works well but labelling features has proved harder than expected.


texts: Displaying geometry points from database layer using QGIS?

texts:
I have an Access Database with a table of Lats and Longs (as well as other attributes).

I can create a layer using an ODBC connection and can open the attribute table etc but can't figure out how to display the points on the screen.

I could create a shapefile of the points separately and link this to the Access file but would prefer not to have to do this.

Is there an easy way of doing this?


texts: How to open .shp file for editing some latitiude and longtitude co-ordinates?

texts:
I have a shapefile. I want to just check data co-ordinates latitude and longitude values in .shp file.
I have using QGIS -2.8.6 , and .shp file automatically loaded in QGIS and it will show polygon output for given co-ordinates.

But, I want to see what is in .shp file. I want the polygon co-ordinates which is shown as a polygon in QGIS. Can you tell me how to see my polygon co-ordinates if the format is in .shp? How can I extract those co-ordinates in CSV or any format?


texts:
I have the following problem:
I have a list of file paths. The files in these paths need to be added to ArcMap. I would like to add these via Python. I know that ArcGIS Pro has a command addDataFromPath() so I loop over the file paths and add them in the command, like so:
for path in pathlist:
    m.addDataFromPath(path)

I've searched around, but I cant seem the find the same command for ArcMap.
Is this even possible? Or is there some other way to load in the files automatically other than simply doing it by hand?


texts: ArcMap runs custom script once and then it crashes on second run

texts:
So I made a script that takes data from an MSSQL database and inserts it into a featureclass. I added the script to a toolbox and am able to run the tool from within ArcMap 10.1.  The problem is that I cannot re-run the tool.  When I try and run it a second time, ArcMap crashes.  I'll start ArcMap again, run the tool and when I try and run it the second time ArcMap crashes again.

By calling arcpy.AddMessage() coupled with time.sleep(5) I have narrowed the problem down to one line of code:
pyodbc_cursor.execute("{call dbo.GetData(?)}",(DataID))

It seems like even though I close() and del all my pyodbc variables (connection, cursor, row) after I'm done with them the cursor is causing a problem when the script is executed the second time around.

This seems like some kind of bug between ArcMap, arcpy and pyodbc


texts: Convert LR/LOLA data to latitude and longitude

texts:
I have height map data from here LR LOLA DEM from the Moon which seems to be "simple spherical" and I want to convert them to latitude and longitude with this:
(It's a bit pseudo code)

from affine import Affine
import rasterio
from pyproj import Proj, transform, Transformer, CRS

infile = rasterio.open(GEOTIFF)
p = (0, 0)  # pixel position

srcProj = Proj(infile.crs)
destProj = ???  # what goes here?

transformer = Transformer.from_crs(srcProj, destProj)

p = p * infile.transform

lo, la = transformer.transform(p[0], p[1])

What destination projection should I use? I tried different ones, but got various exceptions (for example: src and dest are not the same celestial body, etc.)

Addendum: Python describes the infield.crs with the following, maybe that helps:

PROJCS["SimpleCylindrical Moon",
GEOGCS["GCS_Moon",DATUM["D_Moon",SPHEROID["Moon",1737400,0]],
PRIMEM["Reference_Meridian",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Equirectangular"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["standard_parallel_1",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]


texts: Setting Attribute form window resizable in QGIS 3?

texts:
I have created/designed an attribute form window in the Layer Properties - Attribute Form using "Drag and drop"-designer for a polygon shapefile.
UPDATE: If I use "Autogenerate", the Attribute form window is resizable.
However, in that case, I can not benefit from several columns leaving me with a lot of scrolling. I have compared the 2 different qml-files (one saved with "Drag and drop" and one saved with "Autogenerate"). But I'm not able to detect where/how/if it is possible to tell QGIS that the window should be resizable within the qml-file.
The number of fields (attributes) will differ for different shapefiles. It may be 10, it may be 50.
I have organized the attributes by collecting the Species into a Group Box called "Arter" with 4 columns. (It looks a little bit strange since the "default" qml-file was made with ALL species, but in most cases, 8-50 species occurs, and the columns are populated from the left and then downwards in QGIS)

Since the number of species present varies, the Attribute form window that appears when digitizing can get quite big, actually too big even if I set to 5 columns in the Group box "Arter". The Cancel and the OK button is not visible on the screen and I'm unable to resize the window.

Is it possible to:

tell QGIS the attribute window should be resizable
set the width of the Text Edit widget and the Range widget to reduce
the width of the "columns" in the Group Box "Arter"

in the Attribute Form or in the qml-file?
E.g. somewhere near/in the qml-file:

    <attributeEditorContainer name="Arter" groupBox="1" visibilityExpressionEnabled="0" visibilityExpression="" showLabel="1" columnCount="5">

Or setting the width for TextEdit and Range?
The range is from 0-7 and the widget could be more narrow: one digit doesn't take much Space.
<field name="nye_arter">
  <editWidget type="TextEdit">
    <config>
      <Option type="Map">
        <Option value="false" name="IsMultiline" type="bool"/>
        <Option value="false" name="UseHtml" type="bool"/>
      </Option>
    </config>
  </editWidget>
</field>
<field name="Akeleie">
  <editWidget type="Range">
    <config>
      <Option type="Map">
        <Option value="true" name="AllowNull" type="bool"/>
        <Option value="7" name="Max" type="int"/>
        <Option value="1" name="Min" type="int"/>
        <Option value="0" name="Precision" type="int"/>
        <Option value="1" name="Step" type="int"/>
        <Option value="SpinBox" name="Style" type="QString"/>
      </Option>
    </config>
  </editWidget>


texts: Mapnik TextSymbolizer showing nothing above OpenStreetMap base layer

texts:
I'm trying to use Python to create a TextSymbolizer to display a trail
name above the trail LineSymbolizer. Below this layer is the OSM base
layer which contains its own data for that trail.

I can display the line but I cannot display any text. I've run out of
ideas, Google is proving fruitless so I'm looking for new ideas, if I
may be so bold.

I've compiled Mapnik with debug turned on only to find that there's no
logging near the text symbolizer code.
I had a go with gdb but I haven't spent enough time to workout how to
debug from python into the Mapnik shared library. If anybody has any
guidance for how to do this then I'd be grateful. I'm happy enough to
debug but I'm not familiar with gdb.

I also edited the osm.xml and took away the trails and trail names
which might be interfering with the TextSymbolizer but I didn't see
any difference.

Relevant code:

    imgx = 400
    imgy = 400

    m = mapnik.Map(imgx,imgy)
    m.buffer_size = 1000

    workoutParams = dict(host='localhost',port='5432',user='postgres',password='password',dbname='postgis')
    workoutParams['estimate_extent'] = False
    workoutParams['extent']='-180,-89.99,180,89.99'
    workoutTable = "(select id, geometry, 'really long trail name' as trail_name from pub_trail_segment where pub_trail_id = %s) as route" % idNumber
    workoutParams['table']= workoutTable
    workoutDatasource = PostGIS(**workoutParams)
    workoutLayer = Layer('route','+proj=latlong +datum=WGS84')
    workoutLayer.datasource = workoutDatasource
    workoutStyle = mapnik.Style()
    segmentStyleRule = mapnik.Rule()
    segmentStyleRule.symbols.append(mapnik.LineSymbolizer(mapnik.Color('blue'),1)) #50 for print
    segmentTextSymbolizer = mapnik.TextSymbolizer('trail_name', 'book-fonts', 20, mapnik.Color('blue'))
    segmentTextSymbolizer.halo_radius = 1
    segmentTextSymbolizer.label_placement = mapnik.label_placement.LINE_PLACEMENT
    segmentTextSymbolizer.allow_overlap = True
    segmentTextSymbolizer.avoid_edges = False

    segmentStyleRule.symbols.append(segmentTextSymbolizer)
    workoutStyle.rules.append(segmentStyleRule)
    m.append_style('route',workoutStyle)
    workoutLayer.styles.append('route')
    m.layers.append(workoutLayer)
    m.zoom_to_box(bbox)
    im = mapnik.Image(imgx,imgy)
    mapnik.render(m, im)

What am I doing wrong?

Many thanks for any ideas

Cros-posted to StackOverflow


texts:
I'm writing a plugin that extracts some data from a Postgres database, saves a point shapefile of the data and displays it in QGIS. As part of that process, I want to generate a QML style file.

I have this basically working - with the exception of the label placement. I want the label placed NE (above and right) of the point, offset X=1, Y=0.

I have the following code. It all works - but the label is placed 'Around point' rather than offset up and to the right.

    #Add labels.
    layer_settings = QgsPalLayerSettings()
    layer_settings.fieldName = "gr"
    layer_settings.isExpression = True

    layer_settings.OffsetType = QgsPalLayerSettings.FromPoint
    layer_settings.OffsetQuad = 'NE'
    layer_settings.OffsetXY = '1,0'
    layer_settings.OffsetUnits = QgsUnitTypes.RenderUnit.RenderMillimeters
    layer_settings.enabled = True

    text_format = QgsTextFormat()
    text_format.setFont(QFont("Arial Black", 6))
    text_format.setSize(6)
    layer_settings.setFormat(text_format)

    labeling = QgsVectorLayerSimpleLabeling(layer_settings)
    layer.setLabelsEnabled(True)
    layer.setLabeling(labeling)
    layer.triggerRepaint()

    # Save the style file for future use in QGIS.
    outputFile = os.path.join(self.SearchOutputFolder, searchCategoryFileName) + '.qml'
    layer.saveNamedStyle(outputFile)

I've had to guess how to use OffsetQuad and OffsetXY. I just can't seem to find any documentation on how to call/set them.


texts: Creating catalog of ModelBuilder models?

texts:
I have to create a catalog with several models that was made with ModelBuilder.

Do you know a way to record basic information of a model and do it several times. i.e. (Path, Variables, elapsed time)?

Something similar to the report generated in ModelBuilder but without the necessity of opening the model because there are a lot of models.


texts: What are calculation methods to remove moire on raster grid?

texts:
Here's a raster that I calculate in PostGIS.

There are points of residential bulidings with population values P, and businesses with some value to them, say X. I calculate

  Sum(X) / Sum(P) for points within 500 m

And the map has artefacts that distract readers. They ask what that is and distrust the map.

(The red figures show areas that are uninhabbited and have no businesses there, but together the quotent increases, because the denominator drops.)

The problem is that whatever way I smooth it, a quotient has same problem with artefacts. Like if I smooth them by distance:

  Sum(X/(500-distance)*500) / Sum(P/(500-distance)*500) for points within 500 m

Of course, the distance factors in the numerator and in the denominator cancel each other.

Are there ways to fix this?

So far, I chose calculating averages by quarters:


texts: Rasterize polygons using gdal/ogr python API

texts:
My understanding so far is that 'gdal' is used for raster data and 'ogr' for vector data. Given this I have been trying to rasterize a polygon data set (output raster metadata is given by a raster template). I have tried with the following piece of code, but I am afraid that there are quite a few things that are probably wrong.
from osgeo import gdal
from osgeo import ogr

#Geodatabase and Polygon Feature Class name
path_gdb = r'Z:\......\Flood_Defences.gdb'
fc = r'FloodDefences'

#Reading polygons with ogr
driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(path_gdb, 0)
flood_def = gdb.GetLayer(fc)
sr = flood_def.GetSpatialRef()

#Path to the tif file where the rasterized polygons will be stored
out_path = r'Z:\.....\raster.tif'

#Template raster and its information
template_path = r'Z:\.......template.tif'
raster = gdal.Open(template_path)
gt = raster.GetGeoTransform()
ndv = raster.GetRasterBand(1).GetNoDataValue()
proj = raster.GetProjection()
driver = raster.GetDriver().ShortName
cols = raster.RasterXSize
rows = raster.RasterYSize
count = raster.RasterCount
dtype = gdal.GetDataTypeName(raster.GetRasterBand(1).DataType)

#Rasterize polygons
gdal.Rasterize(out_path, gdb, format=driver, outputType=gdal.gdalconst.GDT_Float32,
               creationOptions=['COMPRESS=LZW', 'BIGTIFF=YES', 'blockxsize=256', 'blockysize=256'],
               outputSRS='EPSG:29902', layers=fc, width=cols, height=rows, xRes=gt[1], yRes=gt[5],
               noData=ndv, targetAlignedPixels=True, allTouched='YES', burnValues=1)

I get the following error:
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-177-7fcf99466a93> in <module>
      3                outputSRS='EPSG:29902', layers=fc,
      4                width=cols, height=rows, xRes=gt[1], yRes=gt[5], noData=ndv,
----> 5                targetAlignedPixels=True, allTouched='YES', burnValues=1)

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in Rasterize(destNameOrDestDS, srcDS, **kwargs)
   1137
   1138     if _is_str_or_unicode(destNameOrDestDS):
-> 1139         return wrapper_GDALRasterizeDestName(destNameOrDestDS, srcDS, opts, callback, callback_data)
   1140     else:
   1141         return wrapper_GDALRasterizeDestDS(destNameOrDestDS, srcDS, opts, callback, callback_data)

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in wrapper_GDALRasterizeDestName(*args)
   3488 def wrapper_GDALRasterizeDestName(*args):
   3489     """wrapper_GDALRasterizeDestName(char const * dest, Dataset srcDS, GDALRasterizeOptions options, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset"""
-> 3490     return _gdal.wrapper_GDALRasterizeDestName(*args)
   3491 class GDALBuildVRTOptions(_object):
   3492     """Proxy of C++ GDALBuildVRTOptions class."""

TypeError: in method 'wrapper_GDALRasterizeDestName', argument 2 of type 'GDALDatasetShadow *'

Also, neither sr nor proj seems to give the information I need about the reference system. In the end I had to hard-code it in the function call ('EPSG:29902').


texts: Disaggregating aggregate GIS data using remote sensing

texts:
I haven't really done much with remote sensing imagery. I have census data for an area, and a TIFF raster image of the area.

I'm sure this is possible, is there a simple way to disaggregate the census data to a "per/household" data?

so if let's say the census has a population of 1,000, distribute the population based on the house footprint for example, or lot size


texts: Saved layouts disappeared from composer in QGIS 3.2.1

texts:
I made editions in the projects created in QGIS 2.1.8 using QGIS 3.2.1 and then edited it again in QGIS 2.1.8. After this I realized that all saved layouts disappeared from composers in all projects that I edited in both versions of QGIS.
Does anyone know if it is possible to restore disappeared layouts?


texts: How to set up Visual Studio Development Environment for ArcGIS 10.3.x to develop Add-In?

texts:
How should I set up a development environment for making add-ins for ArcGIS 10.3? There is no walk-thru anywhere that I can find for this. First, ESRI says Visual Studio 2013 is supported but in other places people say it is not. Also here they walk-thru fixes for the MSBuild targets which is a problem I am encountering because I am updating an Add-In but can't even find the targets to rename after installing the SDK.

Furthermore, where is the best link to the SDK to download for 10.3.x development? There is no help here with with what version of the SDK to download as I would expect. What is the order in which the SDK and the best version of Microsoft VS needs to be installed in? I reverted to VS 2010 and removed and re-installed the SDK ArcGIS Runtime for .NET 10.2.6 from the developers.arcgis.com site and there are no targets in my Program Files (x86)/MSbuild/ directory as apparently there should be.

In summary, what is a step-by-step walkthru for developing Add-Ins for ArcGIS 10.3+ with .NET and Visual Studio (any version)?


texts: How to make scripts visible for import in PyQGIS

texts:
This script imports the desired module when used outside of a PyQGIS  script in QGIS.
How to tell QGIS where to look for the scripts, if not with sys.path.append()?

import sys
import os

try :
    import feature_switcher
except:
    print("raw import failed")

print(sys.path)

path_to_scripts = os.path.join('C:\\','Workplace','row_detector','src','qgis')
print(path_to_scripts)
sys.path.append(path_to_scripts)
try :
    import feature_switcher
except:
    print("import with added Windows_style sys.path also failed")

path_to_scripts = 'C:/Workplace/row_detector/src/qgis/'
print(path_to_scripts)

sys.path.append(path_to_scripts)

try :
    import feature_switcher
except:
    print("import with added Unix_style sys.path also failed")


texts: Seeking Crime Analysis Solution Tools in QGIS

texts:
I want to create similar maps like shown at the bottom but I am "limited" to QGIS only right now.
My goal is to compare crimes committed in 2019 to those in 2020 in the same date window of quarantine.
ArcGIS Pro has solution for making the analysis of crime data easier called the Crime Analysis Solution.
How can i achieve similar results (maps) using QGIS?
"Using the Density Change Tool that is a part of the Crime Analysis Solution ribbon, I am able to compare the density raster of a current period to that of another period, essentially working as a raster calculator. The output shows the change between the two with areas in red being higher in change and areas in blue being lower in change"
"The Percent Change Tool in the Crime Analysis Solution ribbon creates a spatial join of point features (incidents) to polygon features (neighborhood boundaries, beats, precincts etc.) and compares two periods of time. The output is a layer showing areas of higher change in red and negative change in blue."


texts: How to label lines with their length in GeoServer?

texts:
I would like to display land parcel dimensions on a map layer. The source layer in GeoServer has the parcel boundaries as line features.

I made an attempt using the geomLength function in an SLD, however the displayed labels appear to represent the length in pixels, not metres or any other real-world unit. The values double when I zoom in and halve when I zoom out.

This same issue was reported a few years ago in this GeoServer forum post, but there appears to be no response.

The label code I'm using is here:

<ogc:Label>
  <ogc:Function name="geomLength">
    <ogc:PropertyName>SHAPE</ogc:PropertyName>
  </ogc:Function>
</ogc:Label>

You can see a sample output WMS map image here

Is there a way to label the lines with lengths in real-world units?


texts: Georeference problem including EPSG 102067 projection

texts:
I'm Qgis user from Czech Republic.
We have a spatial national coordination system (Krovak system) having positive X axes pointing to the south and positive Y axes to the west.
Many geographical systems have problem to display such coordination system, so usually is used transformation which gives:

X_transformed = - Y_Krovak
Y_transormed = - X_Krovak

The above described transformed coordination system is described in EPSG database as EPSG: 102067

To my sadness this projection was excluded from Qgis projection database starting Qgis 1.7.4. Many Czech wms servers  work with this projection and I need their data. Many GIS projects are accomodated into this transformed coordination system.

To bring EPSG 102067 projection into my QGis 1.8 I override its original projection database file which is distributed with projection database file of Qgis 1.7.3 (including EPSG 102067 projection). Name of this file is srs.db

My problem is that I can not correctly geoference raster map into EPSG 102067 projection system. Above mentioned information possible could be part of the problem, or absolutely not. My attempts are resulting into georefferenced raster having height and width less than 1 meter.

Georeference points and original bmp raster could be downloaded form http://www.1progeo.cz/download.php

Login: qgis
Password: qgis

as a zip archive so to debug my mistake could be quite easy. Would you be so kind to help me to give me guidance how to georeference raster map correctly?


texts: How to improve ArcGIS Desktop 10.2.2 and Postgres 9.2 PostGIS geometry performance?

texts:
I'm experiencing sever performance issues on a polygon featureclass(PG_GEOMETRY) in a PostgreSQL 9.2 database from a 10.2.2 client.  Even when zoomed to an area with  10-20 features within the extent.

I've narrowed it down to a call from ArcGIS Desktop:
DECLARE sdecur_508_23777 BINARY CURSOR WITH HOLD FOR  select  st_asewkb(ST_setSRID(zzzztablenamezzzz.shape,-1)) AS shape  from  sde.fooschema.zzzztablenamezzz
and then the  subsequent fetches (e.g. 2015-03-13 19:52:26 GMT LOG statement: FETCH FORWARD 1000 from sdecur_508_23777)

The data frame and the spatial reference id (54017) are the same, so I'm not seeing why a query that brings back the entire table's geometry as SRID-less WKB  needs to be executed (count appox 900000 polygons).

There's no fancy rendering--I just added the featureclass and zoomed in to a small extent.  The featureclass has a spatial index that performs well


texts: Geoserver: How to enable REST Authentication filter chain to allow PUT/POST calls

texts:
As a followup of this question I am wondering how to enable a certain user to be able to place HTTP PUT/POST calls to /geoserver/rest/styles/ without granting the user admin rights or disabling chain security check under Authentication->Filter chain.

With my $geoserver_home/data/security/rest.properties settings being

/rest/**;GET=ROLE_REST_READ
/rest/**;POST,PUT,DELETE=ROLE_REST_WRITE

I assigned user 'chief' the roles 'ROLE_REST_READ' and 'ROLE_REST_WRITE'.

When doing PUT/POST calls with this credentials it will result in Status Code 405 NOT ALLOWED. GET calls however work fine. Using admin credentials or disabling REST chain security check PUT/POST calls work fine (Status Code 200). I tried changing the REST chain settings (check HTTP method, role filter) with no luck. Changing the admin role of the current role service 'default' from ROLE_ADMIN to ROLE_REST_WRITE also works, giving user chief complete admin right over the web interface however.

How/where can I edit the REST filter chain settings to not block PUT/POST calls? Please excuse me if I am missing some basic filter chain points here as I am pretty new to this topic.


texts: How to create LineString geometry in QGIS3?

texts:
The default geometry type for a newly created line seems to be MultiLineString.  It can be changed to LineString using the Multi to Single tool while saving in memory.  However, once it is exported it saves as MultiLineString.  Is there a way to create a line .shp as LineString.

QGIS 3.2.2 Bonn
Mac os X 10.13.6


texts:
I am working on a really big script right now where I have a csv file that I have removed rows and columns from, and edited the headers. I need to create one big shapefile for the entire csv file then create individual shape files for the units under one of the headers. I thougt the best way to do this would be to use arcpy.MakeXyEventLayer(), I saw in an arcgis sample script to then use arcpy.GetCount()for the output file of the xyEveveLayer, then arcpy.SaveToLayerFile_management() and arcpy.FeatureClassToShapefile_ conversion, but when I run the script only my csv file is getting edited and there is no layer in the output file. Is there a step I am missing or should this be making my shape.

this is the few lines of code I have used, after all of the csv file editing, to do what is described above:

outLyr = sys.arg[3] # shapefile layer output name
    XYLyr.newLyr(csvOut, lyrOutFile, spRef, sys.argv[4], sys.argv[5]) # x coordinate column; y coordinate column
    print arcpy.GetCount_management(lyrOutFile)
    csv2LYR.saveLYR(lyrOutFile, curDir)


texts:
I have made a List containing all the original data (.shp files), and I have another List containing the modified data. I basically want to check all the modified file data against the original... So here lies the problem, I next a loop within a loop to check if the original exists within the modified List... and it always seems to find a match? Is my approach right?

This is the sample code I have for PointZ type shapes

 foreach (PointZ origrecord in origPointZList)
        {
            Boolean foundIt = false;

            foreach (PointZ modrecord in origPointZList)
            {
                if ((origrecord.X == modrecord.X)
                    && (origrecord.Y == modrecord.Y)
                    && (origrecord.Z == modrecord.Z)
                    && (origrecord.M == modrecord.M))
                {
                    foundIt = true;
                    break;

                }
                else
                {
                    foundIt = false;
                }

            }
            if (!foundIt)
            {
                result += "Record not found";
            }


texts: Using QGIS3 Processing algorithms from standalone PyQGIS scripts (outside of GUI)

texts:
I'm writing a script which must work outside of QGIS GUI. I call some API functions from qgis.core but I'd like to use the processing plugin.

I'm able to import processing with sys.path.append() but I can't run any process. Moreover, all "native" algs are missing in QgsApplication.processingRegistry().algorithms()

So is it possible to run processing that way ? What am I missing ?

import os, sys
from qgis.core import *
QgsApplication.setPrefixPath('/usr', True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
import processing

layer1 = QgsVectorLayer('data/ROUTE_PRIMAIRE.SHP')
layer2 = QgsVectorLayer('data/ROUTE_SECONDAIRE.SHP')

processing.run('qgis:union', layer1, layer2, 'test.shp') # returns nothing

I'm on QGIS 3.0.1 - Debian 9


texts: Correcting offset between LiDAR datasets

texts:
I am analyzing multiple LiDAR datasets looking at the changes between pre-wildfire and post-wildfire topography. The datatsets are of the same general location. I did not acquire the datasets on my own, so I am not the one who originally reduced the raw data from the instrument, which is part of my issue as I was not the one to verify the horizontal control.

The problem that I am having is that the datasets are, in general, slightly offset from one another. This poses a major problem for generating a comparison (raster subtraction) looking at elevation change. When I process such a comparison, the offset is more apparent on south facing slopes because they are steeper. I have included an illustration of my issue below. As is evident, the south facing slopes show a characteristic false positive change in elevation and the north facing slopes show a characteristic false negative change in elevation, though it is less drastic because the north facing slopes are not as steep.

I am hopeful that someone with experience can shed some light on a possible solution to this issue. Of course, the datasets can be shifted using the georeferencing toolbar, but I would just be guessing at the correct shift. I am looking for a systematic way of handling the shift.


texts:
I am having issues importing an SRTM HGT file into Grass. I first made a new location as the instructions from this Grass Newsletter Tutorial states. I then went to the download site and downloaded the file to my desktop. Inside GRASS, I went to

File > Import raster data > SRTM HGT import [r.in.srtm]

I choose the filename of the .hgt.zip file and remove the extension .hgt.zip (as the instructions say).

The statement looks like this:

r.in.srtm input=C:\Users\Joe\Desktop\N47E007

The following is the output that I'm getting:

Extracting 'C:\Users\Joe\Desktop\N47E007'...
Converting input file to BIL...
mv: cannot stat `C:\\Users\\Joe\\Desktop\\N47E007.hgt': No
such file or directory
ERROR 4: `C:\Users\Joe\Desktop\N47E007.bil' does not exist
in the file system,
and is not recognised as a supported dataset name.
rmdir: C:\Users\Joe\Documents\grassdata/demolocation/PERMANE
NT/.tmp/8120.0: Directory not empty
ERROR: Unable to import data

At the beginning of the script, the file is definitely there because the program recognizes it and starts converting it. Once it produces an error, I notice the file has been moved to another directory

C:\Users\Joe\Documents\grassdata\demolocation\PERMANENT\.tmp\8120.0

What can I do to fix this problem? Thanks!


texts: Leaflet center Marker AND Popup on Map viewport

texts:
I want center my marker on popup open.. and centering map not in marker latlng, but on center of marker and popup! The problem is that popup has dinamic content(loaded on click).

The map size is full display size in a mobile device! I'm just used autoPanPadding option in popup but not sufficient

Refer to follow picture:


texts: Joining New Column to DBF if OpenOffice is Giving "Longer Than Number of Digits Allowed" Error

texts:
When I need to add a data column to a DBF, I always download the shapefile, open the DBF with OpenOffice for Mac, us a VLOOKUP( ) to add my data to a new column, re-save the file as a DBF, then open in in QGIS and that's never not worked for me.

When I do that now, I'm getting an error. I'm using this TIGER/Line Shapefile of Pennsylvania zip codes (.zip), which works fine if I put it right into QGIS and start working with it -- then I add the new column of data, but the second I try to re-save in OpenOffice, it gives me this error (regardless of whether I add the column or not):

The cell it's referring to is simply "108807," so I can't reduce the decimal points that aren't there.  I tried re-copying and pasting as values and it didn't fix it either.

So my two questions are:

Does anybody know how to fix this so I can continue using my simple/preferred method of OpenOffice?
Or can anybody direct me to another way to join the new column I need onto this DBF if it's currently a csv file that can be matched with a unique identifier on the original DBF? (I read something on here earlier about being able to do this with R, but I couldn't figure it out.)


texts:
My use case is one of defined service area or the inverse of the nearest neighbor query. I have a database of service providers each willing to work within N miles (radius) of a point. On the flip side I have users each at a given point looking for service providers that can provide service to their location. Another way to think about this could be restaurant delivery.
If I specify an address, how do I match it to supported delivery areas?
Some might be have a 5 mile radius, some might have 40 mile radius defined in the DB. I say this is the inverse of the nearest neighbor because each "neighbor" in this case has their own definition of how far they are willing to travel from a point, which is more important than how far each is from the customer's point (address).
I'm looking for some SQL query examples.


texts:
I would like to use an aggregate version of the ST_Intersection function in PostGIS. In a table I have multiple polygons and would like to calculate the intersection of all of them. Ideally it would even be possible to group by one of the other columns.
In short, I would like to use ST_Intersection as an aggregate function, or something similar to

SELECT groupid, ST_Intersection(geometry) FROM tablename GROUP BY groupid;

The result of this would be a table with a line for each groupid. The resulting geometry for each groupid should be the intersection of  all polygons with this groupid.

Is there a way to do this in PostGIS?

Example:
My table contains four records, each containing a polygon. They are shown in the left column in the picture below. In the picture assume that the squares are every time the same square (e.g. [-1,1]x[-1,1]). Each of the records has a groupid, shown in text and with a color. In the output I want a record for each groupid. The polygon in the output should be the intersection of all input polygons from that group. In the picture they are shown in the right column.

In this pictures that means the record for groupid = 1 (green) should contain the intersection of all green input polygons and analogous for groupid = 2 (red). Here, the intersection of both green triangles is a smaller triangle and the intersection of the red rectangles is a smaller red square.

The example in code: The code for this example would look like this.

 CREATE TEMPORARY TABLE simpleexample
(
    id INT,
    groupid INT,
    polygon GEOMETRY
);

INSERT INTO simpleexample VALUES
(1, 1, ST_GeomFromText('POLYGON((-1 1, -1 -1, 1 -1, -1 1))')), -- lowerleft half of square [-1,1]x[-1,1]
(2, 1, ST_GeomFromText('POLYGON((1 1, -1 -1, 1 -1, 1 1))')), -- lowerright half of square [-1,1]x[-1,1]
(3, 2, ST_GeomFromText('POLYGON((-1 0, 1 0, 1 -1, -1 -1, -1 0))')), -- lower half of square [-1,1]x[-1,1]
(4, 2, ST_GeomFromText('POLYGON((0 1, 0 -1, 1 -1, 1 1, 0 1))')); -- left half of square [-1,1]x[-1,1]

-- making the intersections manually works fine
SELECT ST_AsText(ST_Intersection(ST_GeomFromText('POLYGON((-1 1, -1 -1, 1 -1, -1 1))'), ST_GeomFromText('POLYGON((1 1, -1 -1, 1 -1, 1 1))')));
-- correctly returns "POLYGON((1 -1,-1 -1,0 0,1 -1))"
SELECT ST_AsText(ST_Intersection(ST_GeomFromText('POLYGON((-1 0, 1 0, 1 -1, -1 -1, -1 0))'), ST_GeomFromText('POLYGON((0 1, 0 -1, 1 -1, 1 1, 0 1))')));
-- correctly returns "POLYGON((0 0,1 0,1 -1,0 -1,0 0))"

-- doing is all in one go doesn't work, this is what I would like to do
SELECT groupid, ST_Intersection(polygon) FROM simpleexample GROUP BY groupid
-- gives error "function st_intersection(geometry) does not exist"
-- I want this to be the two previously calculated lines in one table.

Obviously, I would have to do this for many more polygons than is handy for copy-pasting.


texts: How to get the coordinate system from a PostGIS database?

texts:
My question is probably amazinlgy stupid, but I can't find the answer myself.

I would like to extract from a database the latitude and longitude of points. For now, I'm getting x(the_geom) and y(the_geom), and I would like to know what is the coordinate system of the results I get.

How can I know this?


texts: I am using leaflet to create a set of points, but I want to view which one was plotted first, last, etc. Is there way to label them such that when the first point is plotted, it will plot a circle with a 1 in it? Then plot the next one with a 2, and continued?


texts: Adding dynamic data from a MySQL server as layer in ArcGIS Online

texts:
I am trying to get my head around ArcGIS Online's Dashboards.
One thing that I am struggling with, is how to add dynamic data that is served from an external MySQL server in CSV format to a dashboard. While it is easy to add the data to the Map Viewer (via "Add > Add Layer from Web"), it appears as if it is not as easy to add the dynamic data source to a Dashboard.
These expect the data to be present as layer that, if accessed via the web, must either be ArcGIS Server web service, KML, WMS (OGC), WFS (OGC), WMTS (OGC). There is also the option document, which simply imports the data as a table but does not make it available to the Dashboard as layer (e.g. to develop a time-series graph).
is it possible to save data from the Map Viewer as layer to be used in Dashboards? Or is there a way to establish layers from dynamically served CSV data?


texts: How to find the total number of people who live in a ward 5 miles or less from a closed bank?

texts:
I am currently using QGIS 1.8.0 and need some assistance with a 'Select by Location' query.

I have a data layer showing all UK wards (containing data on ward names, areas and populations) and another with all bank branches. I have run a query to find all the bank branches that have closed. I now want to run a query to find the total number of people who live in a ward 5 miles or less from a closed bank. Please advise as to how to do this.


texts: How to calculate cell tower polygon and get the coordinates in wkt format

texts:
I have an use-case where I have to calculate the polygon of cell tower and store the polygon information in wkt format. I have got the following information regarding a cell tower.

latitude
longitude
cell radius
antenna direction (azimuth)
radius extension (width in degree)

I am looking for an algorithm which uses the above attributes and derives the polygon, which I can save it in WKT format. The polygon would be of petal/teardrop shape (refer the diagram below). The cell radius forms the length of the petal and radius extension forms the width of the petal in degrees.

It would be helpful, if I can get a pointer to any library in Java or Python which I can use for the same.


texts: Move polygon to new location (find new lat-lng coordinates)

texts:
I need to move my polygon (region of geographical area) with 4 different (O1, O2, O3 and O4) lat-lng coordinates (as shown in the figure) to another place, so that I get new values of all 4 lat-lng coordinates (N1, N2, N3 and N4).

So if I have one lat-lng for example the N1, how can I calculate other three coordinates (N2, N3 and N4)?

I am using Parallelogram formula to find the new coordinates.
For example, if I need to calculate N2 coordinate, I will use Points O1, O2 and N1 as the three coordinates of a parallelogram and use them to calculate N2 coordinate.

But is that a correct way to calculate all coordinates? I feel this could work properly for a XY coordinate system but is it valid for lat-lng coordinates system as well?


texts: Trying to add autocomplete on a search bar linked to a CartoDB map column

texts:
Here's my code for the search bar I added to Carto map, with some added script for my attempt at an autocomplete. So far, it hasn't worked. I wish I could be more precise about the problem, but I really don't know what it is. I'm not a developer and just hacking away from what I've found from other people online.

    function main() {

        var map = L.map('map', {
            zoomControl: true,
            center: [LAT, LONG],
            zoom: 15
        });
        L.tileLayer('http://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png', {
            attribution: '&copy; OpenStreetMap contributors, &copy; CartoDB'
        }).addTo(map);

        cartodb.createLayer(map, 'MY_MAP.json').addTo(map)
            .done(function() {
                $('#searchButton').click(function() {
                    input = $("#ad").val();
                    input = titleCase(input)
                    console.log(input);
                                 function initAutocomplete() {
                                      function log( message ) {
                                        $( "
" ).text( message ).prependTo( "#log" );
                                        $( "#log" ).scrollTop( 0 );
                                      }
                                  }
                    var sql = new cartodb.SQL({
                        user: 'bmm80'
                    });
                    sql.getBounds("SELECT * FROM MY_DATA where DATA_COLUMN ILIKE '" + input +
                        "'").done(function(bounds) {
                        $(".title").fadeTo('slow', 0);
                        map.fitBounds(bounds)
                    });

                });
            });
    }