texts: POSTGIS Dissolve operation not working properly (ST_union)

texts:
I want to dissolve a polygon data set in POSTGIS. So far I've read about some different approaches, e. g.  Is there a dissolve function in PostGIS other than st_union?

When I run one of the following queries:

SELECT (ST_DUMP(ST_UNARYUNION(ST_COLLECT(geom)))).geom from myTable;
SELECT (ST_DUMP(ST_UNION(geom))).geom from myTable;

most, but not all features are dissolved. Here an example how it looks like before dissolving:

And after dissolving:

On the left hand side, dissolving was successful, on the right it wasn't.
Does anybody have an idea why this might happen? I want to dissolve all neighboring polygons.

EDIT: The problem becomes awkward. When I Buffer the Polygons before dissolving, the problematic polygons are still not dissolved:

I really can't explain this. The Polygons stem from a ST_DumpAsPolygons operation on a raster..


texts:
I'm new in QGIS world :)
I would like to have an auto-completion script within the editing form on QGiS.
When editing a point features with QGIS form, I would like the script to identify the list of features (polygons) falling under this coordinates and to put their ID in a dedicated field of the feature which is edited.
I have no idea what would be the easiest way to it without reinventing the wheel!
so the steps are : I start the editing of a point feature in the QGIS Form, the script gets the list of polygon features under the coordinates and fills in a myPolygons fields with their IDs...


texts:
In the picture1 shown there is a linestring which is Area2 and a buffer of linestring which is Area1. I want to include Area2(linestring in picture) in the buffer but I cannot increase the size of the buffer because then I will have other extra street from Area2 as well. I have used ST_Snap as well as ST_SnapToGrid but nothing is working.

Is there some other way in postgis that this weird problem is solved without increasing the size of the buffer?

If I increase the size of buffer I will have small line(Extra Area2 line) in picture2 which I dont want in geometry.

Do I need map matching in this case?

or there is some another technique to get this line inside buffer?

or is there any technique to move both Area1 and 2 to same coordinates?

In the end I want all Area2 lines in the buffer without having the extra lines (extra line means here which is shown in picture2).

Area1 and Area2 geometry shapes are same except there are some extra lines in Area2 which I dont want by comparing with Area1.

PICTURE1 showing the difference between geometries. The blue buffer is Area1 and line is Area2. I want Area2 in buffer

PICTURE2 If I increase the size of buffer the small line is also included which I dont want

Picture3 :Another Example If I increase Buffer size it will have more purple linestring from Area 2 which I dont want


texts: Visual Studio Lightswitch and spatial data

texts:
I'm looking for a way to handle spatial data (for instance geography data type in SQL Server 2008) in Visual Studio Lightswitch 2011.

Preferably with an extension that can show data on a map. So far the only extension I have found is the one from NetAdvantage.

Is there alternatives to this extension? Are there (open source) projects addressing this problem?


texts:
QGIS Server seems to be fairly simple to setup and create services, however I'm wondering if there are any limitations with using it opposed to MapServer or GeoServer.  My current server configurations is:
System - One GIS Server

Environment - GIS, DB Server on same machine
Users - 50-100 internal users per day
Function - serving data via OpenLayers for data viewing and data extraction
Using QGIS Server are there any obvious limitations or shortcomings with the following:

User manual/general help
Creating WMS/WFS/WFS-T
General Admin. tool/interface of the server for setting roles, security..etc
Creating cached WMS tiles
Serving out large PostGIS tables (million plus records per table)
Updating services


texts:
I'm trying to make a Flex application where there's going to be 2 main users [Administrator,User] The Administrator can add figures to the map, like polygons, or circles to remark the ratio of some locations, and also notations.
Ex. Mappoint in California, other in Washington, Other in New York etc. Then the Administrator can draw circles in the locations and add notations like " secure zone X km"
And I want the user can see that map with that last modification but he cannot modify, only watch the map the circles and notes that the admin added previously.

I already code something with the code examples from the FLEX API page, but I can't find an example of making a map, configure it, save it, publish and get it later.

Do I need credits from ArcGIS?

Does anybody have a an example, link or something more helpful than drowning in the internet ?


texts: Count gaps between polygon in spatialite

texts:
I have a database (.db) with a table that contains thousands of polygons containing id and geometry. I like a query that can get me the 'number of errors where there are gaps between polygons".

How can I prepare such a query ? What GIS function I can use ?

let's suppose table structure ::

table name : temp1

column name ::
id (integer) : POLYGON_ID
geometry


texts: How to represent the time(weekly/monthly) variationof a polygon map using one image?

texts:
I have a base map made up of three hundred polygons,

Each polygon has lots of attribute values(For example, each week/month get a value). But I only want to draw one picture(due to page limitation) to represent the variation info.

I have tested two methods, one is layers, one is "subplot"(MATLAB).

1.Every monthly attribute can be drawn as a layer, so the picture contains 12 layers. But the result is not clear for some areas, some polygons are covered by others.

2.Draw subplot, a page can contain 12 subplot images, but 52 (week number of a year) is too much for one picture.

I could use both two method with my monthly value, but none for weekly value.
So my question is:
 Is there any other method for visual representation.

I can't paste my own data due to secrecy agreement, but here is example in ArcScene using ArcGIS sample data.( Only three layers was used in sample case, but my own data has 12 or 52 layers. I want to show all the layers using less even only one picture.)

PS: I don't know if there is any better solution except above two method.
Need your help.


texts: GeoServer: Cross layer filtering with current bbox

texts:
I have two layers "building" and "soil" on a French region and I have about 1 million land and 500,000 buildings. I am trying to filter on the WMS query to get the land with the individual type buildings. I added in my CQL_filter

INTERSECTS(wkb_geometry, collectGeometries(queryCollection('poc:building', 'wkb_geometry', 'type=''individual''')))

Unfortunately the number of data is too large and therefore it exceeds my memory limit.
Is there a way to add the current BBOX so that it only filters on displayed entities?


texts: Changing order of fields on CSV export using ArcPy with ArcGIS Pro

texts:
I'm exporting a feature class as a CSV through Python with ArcGIS Pro. I've listed out all the fields I need to keep on export.
Is it possible to use this as a list for ordering my headers in the CSV?
There's a csv.DictWriter but I can't figure out where/how to combine it. I thought listing it out preserves the order in export but it does not.
fields = arcpy.ListFields(fc3)
fieldsToKeep =['REPORTER','RECOWNER','SCINAME','COMNAME','RECBASIS','GlobalID','OCCSTATUS','MANAGESTAT',
               'POPSTAT','OBSDATE','COUNTRY','STATE','COUNTY','COMMENTS']

field_names = [field.aliasName for field in fields if field.name in fieldsToKeep]

with open(outputCSV,'w', newline='') as f:
    w = csv.writer(f)
    #w = csv.DictWriter(f, field_names)
    w.writerow(field_names)
    for row in arcpy.SearchCursor(fc3):
        field_vals = [row.getValue(field.name) for field in fields if field.name in fieldsToKeep]
        w.writerow(field_vals)
        del row


texts: Drawing more complex graphical objects in OpenLayers?

texts:
I'm beginner in OpenLayers and I'm wondering how to draw more complex graphical objects such as the following:

I've taken a look at their examples, but not found that kind of drawings. I'd like to use their primitives, not images.

I'm able to draw the line and put the text on top of it, but I'm not sure how to add small lines at the start and end points. Please, pay attention these small lines are drawn with some angle (not 90 degrees).

The code I've written so far looks as following:

drawSize (size) {
  const draw = new Draw({
    source: this.vectorLayer.getSource(),
    type: 'LineString'
  })
  this.sizeDraw = draw
  this.devDocMap.addInteraction(this.sizeDraw)
  draw.on('drawend', this.onDrawEnd)
},
onDrawEnd (event) {
  this.sizeFeature = event.feature
},
onSizeSet (title) {
  this.sizeFeature.set('graphics', true)
  this.sizeFeature.set('style', sizeStyleId)
  this.sizeFeature.set('title', title)
  this.setFeatureStyle(this.sizeFeature)
},

setFeatureStyle (feature) {
  const styleId = feature.get('style')
  let style = null
  switch (styleId) {
    case 0: {
      style = this.getRedPoint()
      break
    }
    ...
    case 11: {
      const title = feature.get('title')
      style = this.getSizeStyle(feature, title)
      break
    }
  }
  feature.setStyle(style)
},

getSizeStyle (feature, title) {
  if (title) {
    return new Style({
      stroke: new Stroke({ color: 'black', width: 1 }),
      text: new Text({
        text: title.toString(),
        font: '18px "times new roman"',
        overflow: true,
        placement: 'line',
        textBaseline: 'bottom',
        fill: new Fill({ color: 'black' }),
        stroke: new Stroke({ color: 'black', width: 1 })
      })
    })
  }
},

I suppose I should keep writing code in getSizeStyle function, but I'm not sure how to add those 2 small lines.

As you can see I'm drawing the main line using interaction. When I make double click on end point, the process of drawing is finished and 2 small lines should be added automatically.


texts: Copy a layer to a different SRS using GDAL

texts:
I'm trying to create one shapefile from another existing one. However, I need to have the SRS of the new shapefile different from the existing. The code is in C# and looks as:

var layer = inputShapefile.GetLayerByName(_layerName);
var option = String.Format("DST_SRSWKT=EPSG:{0}", WebMercatorWkt);
var options = new[] { option };
outputShapefile.DataSource.CopyLayer(layer, "ShapefileName", options);

The code, however, crashes on the last line without throwing any exception etc whatsoever.

Can someone explain how can I use the CopyLayer method to create a new layer in a new shapefile?

Is there any documentation of what options can be provided as part of CopyLayer() method when copying a layer in the shapefile?


texts: Leaflet: click through marker's transparent background

texts:
The icon I use on my marker does not cover the whole quadratic area.

Is there a way to be able to click through the transparent parts of the marker onto another marker or the map itself?

I'm using only SVGs as marker icons. So I thought of css pointer-events, but this didn't seem to work.

https://developer.mozilla.org/en-US/docs/Web/CSS/pointer-events


texts:
I'm trying to generalize my line features with QGIS. I've tried the "Generalizer"-Plugin, but can't get the results i want. I'd like to generalize my line as shown in the image.

I've tried the "vertex reduction", followed by the "Chaiken's Algorithm" (for example input value: level 10, weigth 4). The smoothing is just too weak and not enough generalized. If i add the level to get a higher smoothing, it results in a failure (memory error). Which algorithms or input values could you recommend?


texts: Is it possible to mimic ArcPublisher with PDF or similar?

texts:
I want to provide mapping output in a format that mimics the delivery provided by ArcPublisher and wonder if anybody has any experience in this.

For instance, the PDF format now supports layers which can turned on and off, but can QGIS or GDAL export in such a way to enable this to happen rather than a basic PDF?

Alternatively, perhaps there is another FOSS product which can provide a similar export.  I want to avoid the client having to download and install software if at all possible (I am assuming that they will have a PDF viewer available).

All suggestions welcome!


texts: Using ArcPy to get point from a cursor for input in arcpy.PolyLine()

texts:
I am using ArcGIS Desktop 10.8 Basic. I do not have access to a Standard or Advanced license for this.
I want to use a point from a search cursor as an input for the arcpy.PolyLine() method. A second cursor will provide the point feature to be used as the and point of the polyline. My search cursor is on a point feature layer, and I'm using 'SHAPE@' as one of the fields in the cursor. This is vexing me to no end. In my mind - and this may be my downfall - there's no reason to drill down to the coordinate level.
Isn't it possible to get the geometry for the current cursor's point feature directly?
How do I actually get from two point features in a cursor as inputs for creating a polyline between them?
OR...does the arcpy.Polyline() method demand an array of coordinates?
I didn't think it did, but it's hard finding good examples. Nearly all the samples and examples I've found in my searches begin with a set of coordinate tuples, never a point feature from a cursor.
Here's the part of my script that is supposed to do the heavy lifting but doesn't work. Some variables are assigned earlier in the script, but I think you'll see my problems don't involve them. Both cursors are on point feature layers.
curfldsAP = ['SHAPE@', 'OID@', fldAPID, fldAPadd]
curfldsGP = ['SHAPE@', 'OID@', "Score", fldGPID]
curNewLine = arcpy.da.InsertCursor(lyrGCline, ['SHAPE@', "FULL_ADDRESS", "ADD_ID"])
with arcpy.da.SearchCursor(lyrAP, curfldsAP) as cur:
    for row in cur:
        linepoints = []
        ptSTART = arcpy.Point(row[0])
        linepoints.append(ptSTART)
        id = row[2]
        add = row[3]
        wc = fldGPID + " = " + str(id)      # whereclause
        with arcpy.da.SearchCursor(lyrGP, curfldsGP, wc) as cur2:
            for row2 in cur2:
                if row2[2] >= 80:
                    ptEND = arcpy.Point(row2[0])
                    linepoints.append(ptEND)
                    newLine = arcpy.Polyline(linepoints)
                    curNewLine.insertRow(newLine, add, id)


texts: How to disable ArcMap Addin combo boxes to not allow extra text next to them?

texts:
I am developing an Addin for ArcMap 10. There are two combo boxes that I use. On several of my users's computers, they get this extra text next to the combo box. It appears it is text from the "Caption" attribute in the Addin xml. I can reset the tool bar but the text comes back after I load any MXD (even new fresh ones with nothing in them). I can not replicate this on my computer, or my co-workers. It only appears on a few users computers in different offices.

I have checked everything they have installed, and found the production tools for 10, but haven't found that causes the problem. They say they have a custom VBA script they use to load some data points, but I can't find any of the VBA code being executed using the VBA editor.

I have also found that in customize mode, I can "Browse for image" and change the image, which makes the text go away, but as soon as I go back to the MXD, this setting is reset. The "Default Style" is checked but greyed out, along with the other options.

Here is my current snipplet of XML for first combo box:

<ComboBox
     id="RSAC_Canopy_Cover.Addin_TargetLayerCombo"
     class="Canopy_Cover.Addin.Combos.TargetLayerCombo"
     message=""
     caption=""
     tip="Target Plot Layer"
     category="Canopy Cover"
     showCaption="false"
     editable="false"
     sizeString="FFFFFFFFFFFFFFFF"
     itemSizeString="FFFFFFFFFFFFFFFF"/>

Here is a picture of the issue:

Has anyone encoutered this problem? Can I disable this in my code, or even overide this in code some how? Does anyone know the cause of this?


texts: NetCDF get_dims too many values to unpack

texts:
I am working with a NetCDF file (.nc) - 600+MB.
import netCDF4
from netCDF4 import num2date
import numpy as np
import os
import pandas as pd

# Open netCDF4 file
file_location = '2m dewpoint temperature.nc'
f = netCDF4.Dataset(file_location)

In order to convert the file to CSV I started by finding its attributes
# Find the attributes
print(f.variables.keys())

Output: dict_keys(['longitude', 'latitude', 'expver', 'time', 'd2m'])

Then, extracted the variable and when trying to get the dimensions
# Extract variable
d2m = f.variables['d2m']

# Get dimensions
time_dim, lat_dim, lon_dim = d2m.get_dims()
time_var = f.variables[time_dim.name]
times = num2date(time_var[:], time_var.units)
latitudes = f.variables[lat_dim.name][:]
longitudes = f.variables[lon_dim.name][:]

I am getting the following error

time_dim, lat_dim, lon_dim = d2m.get_dims()
ValueError: too many values to unpack (expected 3)

What is happening here and how should I solve it?

Edit 1
The output of print(d2m.get_dims()) is
(<class 'netCDF4._netCDF4.Dimension'>: name = 'time', size = 94750, <class 'netCDF4._netCDF4.Dimension'>: name = 'expver', size = 2, <class 'netCDF4._netCDF4.Dimension'>: name = 'latitude', size = 33, <class 'netCDF4._netCDF4.Dimension'>: name = 'longitude', size = 53)


texts:
I just had my mind blown in answer to my previous question. I was going to rework that question, but decided it makes more sense to ask a new one.

I had been converting between ECEF and lat/lon, and ran into an issue with Geoid vs Ellipsoid height. I was looking for an easy way to have an intuitive sense of position on the earth as well as elevation without having to worry about MSL, ellipses and imperfect models.

I would like to define a newLatitude and a newLongitude that are based on a tiny, perfect sphere located at the center of the earth. This sphere rotates with the earth just like ECEF does.

To define a location in this system you must have a triplet {newLat, newLon, Radius}. Using only {newLat, newLon} defines a line radiating from the center of the sphere out into space. You have to say how far along that line you are with the Radius. A radius of 200 meters is 200 meters from the center of the earth. The radius I'm sitting at right now is 6,499,100 meters, etc.

Then we don't have to worry about the average level of the tide for the last 19 years to figure out MSL. We simply use the same reference point for Radius.

Essentially, this is using spherical coordinates for something that only approximates a sphere, but that's fine, because we always have Radius.

To me, this seems to simplify a lot of the problems I have with different coordinate systems. Yes, sure, it leaves unanswered how to define the surface of the earth in this system, but that's where DTED information comes in.

Is this a feasible coordinate system that could be used at all locations on or near the globe without ambiguity? Are there any problems with this scheme (ignoring inertia of existing schemes and reluctance to change... I'm only asking about the validity of the system.)


texts:
I placed some cones on a sports playing field and estimated their locations using latitude/longitude coordinates. An example pair of estimated latitude/longitude coordinates is 54.96975, -1.51407.

On average, my estimated latitude was 0.00009 away from the actual latitude. On average, my estimated longitude was 0.0001 from the actual longitude.

Am I able to say, on average, what distance my estimated coordinates were from the actual coordinates?


texts:
I've imported a SHP via zip folder, and it appears to import correctly in the Data View.

However, in Map View, when I select the Category option one of my twelve items is listed as "Others." In addition, when I change the visualized column, a different item is listed as "Others."

Does anyone have any idea as to what I'm doing incorrectly? Does my question make sense?
I don't know if this has anything to do with it, but I had imported an earlier version, but deleted it. Does it have to do with importing an updated file?


texts: Use user-created GIS layer for online survey

texts:
I'd like to achieve the creation of an online survey where, within a particular study area, different social groups can indicate their use (presence/absence).

In ArcMap (offline) I have a layer containing the study site, which overlays a fishnet with a total number of grid cells of a particular size. I'd like to use these 2 layers as a map in an online survey, where groups can indicate by clicking, per grid cell, if they use the area within that cell.

My questions are:

Whether this is indeed possible and
How I would go on about it and
What survey website can be used for this. I'm relatively inexperienced with GIS and have never done an online survey.


texts: Save a layer for each unique value using Python in QGIS

texts:
I have a map, in which each region is defined as an attribute of the map (I can select each region from the attribute table). I am trying to create new layers each one containing only 1 region using Python (quite like the gui command " save selection as..").

I have got so far (I am very new to Python, so forgive all the obvious errors I have made here):

from qgis.core import *
from qgis.gui import *
layer_Map = canvas.layer(0) #select the layer
provider = layer_Map.dataProvider()
fields = provider.fields()

writer = QgsVectorFileWriter("P:\Qgis\Output\Writer.shp", "CP1250", fields,provider.geometryType(), provider.crs(), "ESRI Shapefile")

Selection = layer_Map.select(1) #This selects the region of the map I want to write in the new layer

writer.addFeature(Selection) #this is clearly wrong, as I get an error


texts: mySQL Spatial ST_contains does not include all Results

texts:
I when I run the first command, mySQL 5.7 only shows a sub-collection of all results. Only when I run the 2nd command I get the correct answer.
The problem with the 2nd command: I cannot use the performance boost of the spatial index anymore (checked with the DESCRIBE command).

SELECT id FROM table1 WHERE ST_CONTAINS(@mls, coords);

SELECT id, ST_CONTAINS(@mls,coords) FROM table1 WHERE ST_CONTAINS(@mls,coords)>0;

coords: Column with x,y POINT spatial data type (with spatial index)
@mls: a user variable containing a multilinestring, made to a polygon with ST_BUFFER.

Does anyone have a solution for that issue?


texts: To select one stream segment at a time for processing

texts: I have attached my stream order network image. I have five different order streams. There are several segments in each order. Now my problem is to select one segment of stream at a time for the further processing. I think this is not a big problem. But I am in trouble, I m using ArcGIS 10 and Python for processing. It would be great if you give some guidance to select each segment of stream at a time using python and GIS10.


texts:
I would like to change a distance unit on Euclidean distance analysis. How can i change the current unit to feet or meter?

what unit is the unit for 0 and 0.18939. is it degree?


texts:
I was reading an article that said when we display a map in WGS84 many people call it no projection. But most definitely when you see the map it has been distorted on a 2D surface.

Which got me to thinking is there such a thing as a "map" with no projection. If I'm looking at a 3D digital globe it is representing how it would look as if I were out in space. So then it wouldn't be projected. But if I'm looking at the globe on a screen then it is being rendered on a 2D screen so it is projected.

Is there such a thing as no projection?


texts:
I want to personalize a form in QGIS 3.4. Just a basic case: when a user open attributes table, a message box fire up. For an incomprehensible reason this message box open a first time, I click "OK" and a new message box appear a second time.

Why this behavior? It's the same with QGIS 3.10.2 (LTR)

from qgis.PyQt.QtWidgets import QMessageBox

def form_open(a, b, c) :
    msgBox = QMessageBox()
    msgBox.setText("Hello")
    msgBox.exec_()

This initialization function needs absolutely 3 arguments, but I don't found any clear documentation about it.

Here are my settings in QGIS:

In layer properties / Form
From a UI file and a Python script
From an external .py file
Initialization function name : form_open

Here is my UI file, with just a QPushButton :

<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
 <class>Dialog</class>
 <widget class="QDialog" name="Dialog">
  <property name="geometry">
   <rect>
    <x>0</x>
    <y>0</y>
    <width>400</width>
    <height>300</height>
   </rect>
  </property>
  <property name="windowTitle">
   <string>Dialog</string>
  </property>
  <widget class="QPushButton" name="pushButton">
   <property name="geometry">
    <rect>
     <x>150</x>
     <y>100</y>
     <width>75</width>
     <height>23</height>
    </rect>
   </property>
   <property name="text">
    <string>PushButton</string>
   </property>
  </widget>
 </widget>
 <resources/>
 <connections/>
</ui>


texts: Determine if Point is in feature class attribute table

texts:
I'm working with ArcObjects, C# and ArcMap. I have an attribute table that has several rows of data including a "shape" field that holds a point geometry. I also have a random point. How can I iterate through each point in the table to see if it matches my random point?

I thought using something like ISpatialFilter would work, but I'm not sure if I'm using it right. My code is below:

IPoint _point = new PointClass();
_point.PutCoords(3605857.580, 516199.563);
ISpatialFilter _filter = new SpatialFilterClass();
_filter.Geometry = _point;
_filter.GeometryField = _featureClass.ShapeFieldName;
_filter.SpatialRel = esriSpatialRelEnum.esriSpatialRelContains;
//_filter.WhereClause = "Shape = "+ _point+"";    // DONT KNOW WHAT TO PUT
//_filter.SubFields = "Shape";

IFeatureCursor _featureCursor = _featureClass.Search(_filter, false);

while ((_featureCursor.NextFeature()) != null)
     {
        if (_point == _featureCursor)
            {
               MessageBox.Show("FOUND");
            }
     }

As you can see, I do have some code edited out as I'm not sure what to put in certain areas.


texts: SQL Server data modelling + spatial data = what book?

texts: What recommended book can you provide about data modelling for geographical spatial data (geometry) in SQL server?


texts: Mapping origins and destinations to the created isochrones in QGIS

texts: I created isochrones in QGIS using Iso-area as polygon algorithm with an interval of 30 minutes. I have some coordinates related to some origins and destinations. I need to map these origins and destinations to the created isochrones to identify what percentages of them lie within that 30 minutes interval. In other words, if I want to travel from the origins to the destinations what percentage of those trips completely lie within the created isochrones? Is there any tool in QGIS to help with achieving this goal?


texts: How does a GPS device determine bearing?

texts:
Following a discussion on Math and this particular comment, I am wondering how a GPS receiver determines its bearing. As I understand it:

GPS is not a directional system, thus bearing can only be inferred from receiver movement. If a receiver has obtained a fix and a valid bearing, then turns while stationary, the bearing will not reflect that until the receiver actually starts moving in that direction.
So far I have therefore assumed that the receiver simply takes the current and previous location, assumes it has traveled in a straight line from one to the other and returns the appropriate bearing.
Bearing will therefore suffer from errors which are linked to the positional accuracy of these two positions, as well as the distance between them.
The error increases whenever the receiver is not traveling along a straight line – rather than the current bearing, the receiver would return the average bearing as it moved from the previous to the current position.

But the comment I got seems to indicate otherwise.

Are my assumptions correct? Or do GPS receivers have another way of determining bearing from satellite signals alone?


texts:
I just wanted to ask around if you have already made any experiences in the usage of the format Geopackage in offline maps (on mobile devices).

Am starting to get into it but can't find much information about the technical implementation of geopackages in map applications. Therefore, would be glad for any hints.


texts: NEAR_3D function without the need for 3D Analyst, Arcpy

texts:
I am looking for a Python script which would do similar as the ArcGIS 3D Analyst NEAR_3D function.

Basically I would like to programmatically in Python find nearest points of one point feature class to another point feature class and calculate their delta X, Y, Z.


texts: Efficient processing of difference between input data and database

texts:
I have an input dataset whose records will be appended to an existing database. Prior to being appended, the data will go through heavy, time-intensive processing. I want to filter out records from the input dataset which already exist in the database to reduce processing time.

The difference between the input and database is illustrated here:

This an overview of the kind of process I am looking at. The input data will eventually feed into the database.

My current solution involves using a Matcher transformer on the combined database and input, then filtering the NotMatched result using a FeatureTypeFilter to retain only the input records.

Is there a more efficient way to obtain the difference features?


texts: How can I create custom projected coordinate system in arcpy?

texts:
I want to create my own projected coordinate system based on one of geographical coordinate systems and set custom values to properties as longitude, scaleFactor etc.

I found SpatialReference class.

But can't understand how set geographical coordinate system that will be base for my coordinate system - GCS property is read only, and no methods for set its value.

And also, when I calls create method, result is RuntimeError: ERROR 999999: Error executing function.

Article has example with only reading existing property, but hasn't examples with creating CS.

What am I doing wrong, and where can I find examples for CS creating?

My code:

new_sr = arcpy.SpatialReference()
new_sr.linearUnitName = u'meters'
new_sr.longitude = 30
new_sr.longitudeOf1st = 10
new_sr.longitudeOf2nd = 50
new_sr.longitudeOfOrigin = 30
new_sr.projectionCode = 100500
new_sr.projectionName = u'my_projection'
new_sr.scaleFactor = 2
new_sr.standardParallel1 = 50
new_sr.standardParallel2 = 60
new_sr.create() # RuntimeError: ERROR 999999: Error executing function

May be I must set some other properties, or vice versa remove unnessecary?


texts: Returning results from script tool in Geoprocessing Results window?

texts:
I have a script tool that prints a list of unique attributes in a feature class. It works fine in PythonWin. When I include it my model, it runs successfully, but does not display the attributes in the geoproccessing results window. I am a very new user to Python and ModelBuilder. I read the help previously on this site "Why are results from Python tool not in results window", but I still don't get it.

Here is my script that managed to put together.

import arcpy
arcpy.env.workspace = "F:\Shared\DDS\GIS\Infrastructure\DATAMIGRATE\CURRENTPROJ\HOUSTON_PLAT_TRACKER\ESRI\ImportAndCleanLines\Scratch\Scratch.gdb"
arcpy.env.overwriteOutput = True
with arcpy.da.SearchCursor ("DwgToUpload", "Layer") as cursor:
    print sorted({str(row[0]) for row in cursor})


texts: I am working on a project to better understand possible applications of the Sentinel 5p dataset within GEE. If anyone has current or past projects that they are willing to share, please feel free to do so. This is out of pure curiosity and academic purposes.


texts: Modifying polygons to be more rectangular using PyQGIS 3

texts:
I'm trying to re-adapt an adaptation of this code but I'm getting the error stated above.
I tried reading the docs, but I honestly didn't understand it.
The input is a vector that contains many polygons and I want to make them more rectangular.
def slope(x1, y1, x2, y2):
        return (y2 - y1) / ((x2 - x1) + 0.00001) #I was getting DivisionByZero error

def intercept(y, slope, x):
        return y - slope * x

def perpindicular(slope):
        return -1 / slope

pth = r"C:\Users\whatever\...myVector.gpkg"
lyr = QgsVectorLayer(pth, "idk what i should put here", "ogr")

epsg = lyr.crs().postgisSrid()
uri = "Polygon?crs=epsg:{}&field=id:integer&index=yes".format(epsg)

rectangle = QgsVectorLayer(uri, 'Rectangles', 'memory')
p = rectangle.dataProvider()

for f in lyr.getFeatures():
    g = f.geometry()
    xmin, ymin, xmax, ymax = g.boundingBox().toRectF().getCoords()
    pts = f.geometry().asPolygon()[0]

for i in range(len(pts)-1):
    if pts[i][1] == ymax and pts[i+1][1] < pts[i][1]:
        idx = i
    if pts[i][1] == ymax and pts[i-1][1] < pts[i][1]:
        idx = i-1
r = []
x1 = pts[idx][0]
y1 = pts[idx][1]
r.append(QgsPointXY(x1,y1))
x2 = pts[idx+1][0]
y2 = pts[idx+1][1]
r.append(QgsPointXY(x2,y2))

s1 = slope(x1, y1, x2, y2)
i1 = intercept(y1, s1, x1)

x3 = pts[idx+2][0]
y3 = pts[idx+2][1]
i2 = intercept(y3, s1, x3)
s3 = perpindicular(s1)
i3 = intercept(y2, s3, x2)
x4 = (i3 - i2)/(s1 - s3)
y4 = s3 * x4 + i3
r.append(QgsPointXY(x4, y4))
s4 = perpindicular(s1)
i4 = intercept(y1, s4, x1)
x5 = (i4 - i2)/(s1 - s4)
y5 = s4 * x5 + i4
r.extend([QgsPointXY(x5, y5),QgsPointXY(x1, y1)])
poly = []
poly.append(r)
g = QgsGeometry.fromPolygonXY(poly)
ft = QgsFeature()
ft.setAttributes([i])
ft.setGeometry(g)
p.addFeatures([ft])

QgsProject.instance().addMapLayer([rectangle, lyr])


texts:
Problem : I want to map using a drone the lake's shores and nearby water shallow areas;
I used dronedeploy to plan the survey ;  I setup the drone to capture at least two pictures overlapping the Shore/shallow water area; The shore area is densely vegetated with trees and shrubs; Shore/water interface is the vegetation or there is a little (few meters) sandy line (2-3m) ;

I fled the drone (mavic pro) to a relative height of 70/75 meters;
Acquisition went perfect and I expected that some of the frames would have been rejected; especially those which present little ununivocal reference mosaic points (ie too much water in the picture)

Unfortunately the result is not as expected and most of the frames (images) are rejected in the stitching process.

Any idea on how to setup a proper drone survey to capture this type of coastal/Shore features?


texts: Is there a recipe for large shapefiles to be used in an Android App?

texts:
I am very new to GIS and I am diving into cartography and mapping because of an Android App that I am currently developing. I want to display environmental data, more precisely noise pollution data on a map in the App. But the data I received is quite complex and large and I am having trouble to properly simplify and export the data that it can easily be used with OpenStreetMap on mobiles.

I have several versions of the data, one older shapefile 23.5 MB, one newer, more accurate one with even 61.3 MB. When I export the files with QGis I get .kml of 66.4 MB and over 100 MB. I failed in simplifying geometry to reduce the file size since the size seems only to get smaller as the shapes stop looking too different from the original.

I heard of the concept to slice/intersect the shape, for example with shapefiles from district borders so the mobile device would only load that part of the large file when being in the respective district, with GPS on. I tried to do this with ArcGIS after having downloaded .kml of district borders, but this software seems too complex for me, unfortunately.

Is this also possibly be done with QGIS?

What is the best way to boil a large layer of data down so it can easily be used for an Android app? I guess these large polygonal shapes can be expressed in a way simpler way to still get the values for a certain location?!

The data, in this case a .kml of 66.4 MB imported via fusiontables onto Google Maps can be seen here: https://www.google.com/fusiontables/DataSource?snapid=S717313eWaJ

I hope some GIS expert has time to look at this and can guide me into the right direction.

Edit 1:
Concerning tiling. I see that tiling works only for raster images. The essence of using the .shp for the is the actual data that lies behind it - the app needs to be able to sort the data, e.g. show only values >x . In order for the device to know on which location which value has been measured would tiling the shape only be a graphical solution?


texts:
I am currently using QGIS 3.14.
I cannot find the OpenLayers Plugin, it is not in the list, and when I add Quick Map Services Plugin, it does not have Bing maps in it (or Google Maps either).
How can I get Bing maps (satellite) to display in QGIS 3.14?


texts: Viewing Multiple Grouped Data as ESRI Map Service in ArcGIS Explorer Online

texts:
I'm trying to add ESRI Map Service to ArcGIS Explorer Online. This map service consists of 3 group layers, one of which contains 3 further subgroups. I tried adding my mapservice to the web viewer using the following URL REST interface structure:

http://myserver/ArcGIS/rest/services/myservicename/MapServer

However, this returns no results. The viewer recognizes the map service, and when I add it up, I can see all the layers and individual shapefiles in the legend window, but not in the actual viewer. It also recognizes the spatial extent for each layer - when I zoom in on a particular layer, the viewer centres on the correct extent. I can turn them on and off in the legend but I am unable to see them.

When I add each individual map service layer separately it works fine (with the exception of a couple of layers), as follows:

http://myserver/ArcGIS/rest/services/myservicename/MapServer/2

The other odd aspect is that the same map service works properly in ArcGIS.com map viewer, even when added by its general map service name (http://myserver/ArcGIS/rest/services/myservicename/MapServer).
Basically, my question is whether you can add an ESRI map service to ArcGIS Explorer Online with multiple grouped data in a single operation, just as in ArcGIS.com?

I'd appreciate any help as I've been unable to find an answer elsewhere.

Thanks,
Magda

P.S. I'm going to post cross post on the http://forums.arcgis.com/


texts: QGIS Modeler different results by calculate the $area

texts:
When I use the graphic modeler,add a field and calculate the area, the sum of the area of all records is different to result which I get by calculation in the fieldcalculator. The result in the fieldcalculator is valid. Alle shp-Files have the EPSG 31468

-->the results by Modeler sum: 10275 m²

-->the results by the fielcalculator sum: 10311 m²


texts: Merge and export two layers that utilize blending modes

texts:
I have spent more time than I can to admit trying to come up with a solution that doesn't require bringing the layers into Photoshop or the GIMP (which would take an insane amount of RAM/swap space to work on this).

Here's the setup:
Using QGIS 2.8.1 on Ubuntu.  2 layers imported, both have srs and align properly.

Layer 1: Low resolution color imagery for the entire planet
Layer 2: Higher resolution B&W imagery for an entire planet

I'm using the color information from Layer 1 and overlaying it on the higher resolution B&W imagery of layer 2.  To do this, I use the following:

Layer one: Properties: Color rendering: Blending mode: Hard Light
Layer two: Properties: Band rendering: Contrast Enhancement: No enhancement
Layer two: Properties: Color rendering: Contrast: -15

After setting these values, the resultant two datasets look very good over one another.

Now, my issues as I see them:

I need to merge the two layers into one
I need to export the resultant merged layer into a geotiff

I can figure out #2 okay using layer: Save As: Output mode: Rendered Image, etc

But, what I cannot figure out is #1.


texts: Saving a gml layer does not change coordinates

texts:
I am trying to export a gml layer to geojson with lat/lng as the preferred coordinate system. So what I did was a 'save as' on the layer with selected CRS: 'IGNF: WGS84G'.

Now I was hoping it converted the coordinates of the GML but it did not.

Part of GML data:

<gml:PolygonPatch>
    <gml:exterior>
        <gml:LinearRing>
            <gml:posList>264124.695 566223.95</gml:posList>
        </gml:LinearRing>
    </gml:exterior>
</gml:PolygonPatch>

Part of converted JSON data:

"geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 264124.695, 566223.95 ] ] ] ] }

As you can see the coordinates remain the same. I tried other formats in order to test if something is converted but not. I am not an expert on GIS transformations so probably I did something wrong.

Environment:
Windows 2012 Standard running inside VMWare Fusion
QGIS version: 2.0.1 - Dafour 64 Bit

A screenshot of my actions:


texts: Problem with custom projection definition in Qgis

texts:
I want to define a custom crs valid for Mexico. The context is that older cartography in Mexico is based on the NAD27 datum. I wanted to test the following proj4 string
+proj=longlat +ellps=clrk66 +towgs84=-12,130,190,0,0,0,0 +no_defs
Nevertheless , when I test my string with 19N 99W (19,-99) I get the following values

testing the proj 4 definition included in Qgis gives the same results.
These are wrong since the resulting north should be around 19 and the resulting east around -99. I wonder if anyone can tell me what I am doing wrong.
Thanks in advance


texts: MapInfo shapefile vs. naturalearthdata.com shapefile

texts:
When I load this map shapefile (exported from MapInfo as a shapefile)

shp.WashingtonLand<-readOGR("C:\\...\\Maps","WLmap2_polyline")

R replies:

OGR data source with driver: ESRI Shapefile
Source: "C:\Users\Andreas\SkyDrive\Speciale\02. Data\02. Statistik", layer: "WLmap2_polyline"
with 2577 features and 23 fields
Feature type: wkbLineString with 3 dimensions
Warning message:
In readOGR("C:\\Users\\Andreas\\SkyDrive\\Speciale\\02. Data\\02. Statistik",  :
  Z-dimension discarded

How do I see what this Z-dimension is?

Can this explain why the map won't show with points on, or as a layer on top of another shapefile like the ones from Naturalearthdata.com?


texts:
I have a map that I want to do graduated chloropleth maps for each year over 10 years, and I would like to use the same graduations for each year.

I know this was asked earlier and therefore I know that I save it as a style and then load it for each new year, but I'm having problems with the save style/load style function in qgis - every time I load the style the data converts back to the original year.

Has anyone else had this problem and figured out a solution? I'm using qgis 2.2.0 valmiera (desktop).


texts: Management of transparency of layer using QGIS Server and mapproxy

texts:
I have problems using mapproxy with QGIS-server via openlayers (mapproxy-1.6.0, qgis-server 2 and openlayers-2.13.1).

I have a layer in which i defined the background of the polygon with transparency.

When I asked directly  qgis-server via openlayers, I get the png wanted (png background is entirely transparent and the transparency of the background is what i need).

When I pass through maproxy, the image background is transparent but not the background of the polygons (the setting of the mapproxy may be ok because I get a transparent png image).
What's weird is when I set the background of a polygon with less than 50% transparency in QGIS, through mapproxy the background of the polygon is completely opaque and when I put a transparency equal to or greater than 50%, the background of polygons is completely transparent.

Exemple directly from qgis-server :

Exemple via mapproxy :


texts: How do I know that ArcMap is running a tool using meters as the unit?

texts: I am not sure how to tell what units ArcMap is defaulting for with the tools it is using. For example, I am currently trying to use the line density tool to create a road density layer. I need to use this layer later on so I need to know the resolution of the output and would like to set it to a specific distance in meters. However, it seems that the tool (and other tools) always default to some other undefined measurement. I have tried changing my document to a projected coordinate system and setting the units to meters so that it actually says "meters" in the bottom right but it seems that the tools are still defaulting to something else.


texts: Create Layer from Selected Features tool: How is the list of features stored?

texts:
In ArcGIS Desktop, there is a tool called Create Layer from Selected Features. Using this tool, I can convert a selection to a layer in ArcMap.

It's my guess that there is a "list of features" that is embedded in the layer somewhere that isolates which features to show in the map.

I would have thought that the "list of features" would have been stored in a where clause in a definition query as a list of objectids. For example: where objectid in(123,456,789,...). However, this is not the case.

If the layer doesn't have a WHERE clause or definition query, then how is the list of features stored?


texts: Does the Clipper plugin for QGIS work with Shapefile layers?

texts:
I created 2 overlapped polygons in one Shapefile layer, run the intersection tool with another layer and then run the Clipper plugin.
This causes the error:

TypeError: MultiPolygon geometry cannot be converted to a polygon.
Only single polygon or curve polygon types are permitted.

Well, I decided not to do intersection before the Clipper. But now I've received the same error...


texts: Copy FGDB from folder to another folder in Python 3 without ArcPy?

texts:
I found this code here. first time I ran it, it copied my file to another folder but for second time nothing happened. Even I didn't get error in Jupyter Notebook.

#Move GDB to G: for your project's folder

src = r'folder_1'
dst = r'folder_2'

import shutil, errno

def copyanything(src, dst):
    try:
        shutil.copytree(src, dst)
    except OSError as exc: # python >2.5
        if exc.errno == errno.ENOTDIR:
            shutil.copy(src, dst)
        else: raise

I don't have access to ArcPy in my PC. I need something which works without ArcPy and for Python 3.


texts: Computer Locking GDB When Running Scheduled Task

texts:
I am currently downloading and unzipping a GDB every night through the use of Task Scheduler. Then I update a related GDB with the changed features in all of the feature classes.

My issue is, the scheduled task script to download and unzip the GDB is creating a lock on the nightly downloaded GDB. This isn't an issue until I run the script on the next night.

My "try and except" statement to delete the previous GDB allows for the script to run but obviously doesn't delete the script do to my computer lock still being on the GDB.

I tried to delete the lock manually through windows explorer but it doesn't work.  It seems that the schedule task manager just holds on to the lock.

This is my first experience with running a script through task scheduler where I actually need to delete an object created by the task scheduler.  the only way I've found to unlock the GDB is by restarting my computer.

Code:
import arcpy, os, re, zipfile, shutil
import datetime as dt
from arcpy import *

#Relevant Dates
today       = dt.date.today()
yesterday   = today - dt.timedelta(days = 1)

#Create Target GDB
print "Creating GDB {0}".format(today)
env.overwriteOutput = True
outLocation         = r'C:\XXX\YYYYY\ZZZZZ GDB'####Update Location with Shared Drive
GDBName             = "Pathfinder {0}.gdb".format(today)
GDB                 = os.path.join (outLocation, GDBName)
arcpy.CreateFileGDB_management (outLocation, GDBName)

#Copied File Location and Unzip
print "Finding Zipped Folder"
sourceFolder = r'C:\XXX\YYYYY\ZZZZZ\Extracts'####Update Location with Pathfinder Data Dump folder
for (dirname, dirs, files) in os.walk(sourceFolder):
    for filename in files:
        if filename.startswith ("Location"):####Update Location with Pathfinder extract
            zipFl = os.path.join(dirname, filename)

print "Copying Zipped File: {0}".format(zipFl)
copyLoc     = r'C:\XXX\YYYYY\ZZZZZ Export'#####Update Location with shared drive copied folder
shutil.copyfile(zipFl, os.path.join(copyLoc, os.path.basename(zipFl)))

print "Unzipping File: {0}".format(zipFl)
zipFl       = os.path.join(copyLoc, os.path.basename(zipFl))
zipOb       = zipfile.ZipFile (zipFl, "r")
zipOb.extractall(GDB)
zipOb.close()

#Clean Up Folders
print "Deleting File: {0}".format(zipFl)
os.remove(zipFl)
yesterName      = "Pathfinder {0}".format(yesterday)
yesterdayGDB    = os.path.join(outLocation, yesterName)
try:
    arcpy.Delete_management(yesterdayGDB)
    print "Deleted {0}".format(yesterdayGDB)
except: print "Pathfinder {0}.gdb does not exist".format(yesterdayGDB)
print "Finished"

ArcGIS 10.1
Intel(R) Xeon(R) CPU E5-1607 0 @ 3.00 GHz
RAM: 16 GB
64 bit OS


texts: I´ve got a problem with displaying the imported gpx file. The gpx file is created in Holux Funtrek 130 PRo. The map I use is from Openstreetmap, I am sure it is the right map from the location I recorded the gpx file. Everything works fine...except the gpx file doesn´t show up. The layer icon is ticked in the layer(s) box. WGS 1984 is set up, both on my gps device and in QGIS. I study geoinformatics, and am really interested in how these things work..but nothing like this has been taught yet..we´ve learnt only adding scalebar and all these elementary tasks. I´d be glad if there was anyone who could answer this question. Thank you.


texts:
I'm trying to print a map using the print manager.  The map I add to the page displays more or less as I'd like it to appear. It is set to zoom to the selected layer, which is a polygon of the area I'm interested in.  It displays in the viewer just fine.

If I go to the print and layout manager, set my page and paper size and add the map, the map content on the page includes a lot of area outside of the selected layer that doesn't appear on the view I'm adding the map from the map.  It does this if I have saved this view as a project as well as if I don't.

I'd like to zoom the content of the page to larger than the page and "crop" it to fit the page or something to display just the parts of the map I want to show.

I can see how to adjust the size of the page to fit the content, but, not the other way around.  The content I get when I add the map is not what it is I think I am adding.

Any suggestions???


texts: Is there any way to add an area column to new shapefile created from raster with gdal's polygonize in the same step?

texts:
I've been going through a fairly lengthy process of polygonizing very intricate raster images, then creating a new layer with the extra area attribute, so that I can then remove the millions of tiny dots that turn what should be a 5 megabyte shapefile into a 40 megabyte one.

When I run the Raster|Conversion|Polygonize with qgis, it shows the gdal script that is run in the command line to complete the task. I don't know many gdal commands other than ogr2ogr stuff, so was wondering if any of you gdal pro's out there knew how to add an area attribute to the outputted shapefile (populated with calculated areas) without a second step by adding something to this script:

gdal_polygonize.py "/Users/douglas/Desktop/polygonized_rasters/zones.tif"

I've clicked a couple of options in this polygonize menu and saw code added to the script, so that gave me the idea that there might be some code I could add to this line to add that area column during the same process. Maybe not, I'm not sure.


texts: Can i use CreateMapSDDraft in the arcpy.mapping library to set the scale range on tiled services at 10.1?

texts: I know that Esri has exposed many properties to the arcpy lib for publishing services. I need to publish a tiled/cached map service but I need to set the min/max scale when I do that. Can this be accomplished via Python code? I don't see it as a property on the help page of the CreateMapSDDraft script. Is it just not supported or am I doing something wrong?


texts:
I have a problem as a short and clear. When I use the snipped code:

 var select = new ol.interaction.Select({
        condition: ol.events.condition.mouseMove
    });

base map's mouse click drag is not working. You can see http://ol3js.org/en/master/examples/select-features.html?q=select

When you select as hover from  Action type. Drag is not working.
is it normal?


texts: I have a vector map with all the zip-codes of Denmark. I have taken a picture in each zip-code, and I now want to fit theses pictures into the zip-code polygons (I'm using QGIS). I can add the pictures one by one by changing the layer style to raster, and scaling the image. Seeing that there are almost 600 zip codes in Denmark, and that the pictures might change through time, I would very much like to find an automated way to do this. Hope somebody can help.


texts: Reading Esri File Geodatabase (FileGDB) with OGR?

texts:
I have a file ending with *.gdb. I don't know what it actually is. I also don't know what to expect inside the file. I am trying to access it with ogr.

I tried with ogrinfo:

ogrinfo test.gdb
ERROR 1: posix_spawnp() failed
FAILURE:
Unable to open datasource `test.gdb' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> GML
  -> GPX
  -> KML
  -> GeoJSON
  -> GMT
  -> GPKG
  -> SQLite
  -> WAsP
  -> PostgreSQL
  -> PCIDSK
  -> OpenFileGDB
  -> XPlane
  -> AVCBin
  -> AVCE00
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDump
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> EDIGEO
  -> GFT
  -> GME
  -> SVG
  -> CouchDB
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> XLS
  -> ODS
  -> XLSX
  -> ElasticSearch
  -> PDF
  -> CartoDB
  -> SXF

So apparently I do have the OpenFileGDB driver. Why do I get an error?

I am on OS X 10.10.3

My GDAL version is 1.11.1, released 2014/09/24 installed with brew.


texts: How to improve legend spacing for ESRI ArcGIS map

texts: I am using ArcGIS 10.1 and am trying to achieve a fairly uniform amount of white space between the first four legend items and am having difficulties pulling this off.  All four items are separate shapefiles/layers and it would be hard for me to turn these into single layers with categories.  I have had success adding space between the two "areas" by changing some of the values in the legend properties, but I cannot seem to figure out how to reduce the space between the two fire station layers. Any suggestions?


texts: Labeling single road once per data driven page?

texts:
I am working with data driven pages and need to label a long road that spans multiple pages.and convert those labels to annotation.

Is there a way for ArcGIS Desktop to treat each page as its own segment, forcing it to place a new label on each page?

I have messed around with the connect features/one label per feature/one label per feature segment/one label per feature part options in the label density tab of the label layer properties but have had no luck.


texts: Reclassify raster values continuously instead of assigning them to specific groups

texts:
I want to create a Habitat Suitability Index (HSI). I have several variables (raster layers containing information about e.g. precipitation, slope, food resources, etc.).
I have to create an "index" for each raster layer. That means that i have to reclassify all layers on a scale from 0 to 1, where 0 is not suitable and 1 is really suitable. This is were I am stuck. I know how to reclassify data in terms of assigning values to specific classes. But for some  variables, I have to create a continuous index
I have troubles explaining, since english is not my first language and I'm  new to the topic, but I will try:
I need to reclassfy the data continuously, according to my guidelines. For example altitude (image 1):
The plot below shows, how i have to reclassify altitude. in summer, an altitude from 1200 to 1500 meter is perfect, for higher altitues the HSI slowly decreases, but not linear unfortunately.
This is the same for example for the index for the crown projection (second image). I have to reclassify the values on a scale from 1 to 100, but again, it is not linear.

My first thought was to apply some kind of funtion on the raster, but unfortunately, I could not find the right answer on google, so I'm afraid I'm searching for the wrong keywords.
I am using QGIS (3.16.1 and 3.10.13), I am quite experienced in SAGA, GDAL, Grass and OTB and I have a bit of knowledge in R.


texts: Automatically updating table column with name of table using MapBasic

texts:
Does anyone know a straightforward way in MapBasic to automatically update a table column with the name of the table the column exists in?

Going from the MapBasic reference guide it seems the TableInfo function should serve this purpose but I can't seem to get it to work. In the MapBasic window whenever you use that function it says not unrecognised command.


texts: Creating Local Cache From ArcGIS Online Basemaps

texts:
I am contemplating creating a local caches from the standard ArcGIS Online Basemaps; streetmap, imagery and topo.  I need to do this because the environment the web application will be hosted is DOES NOT have live connection to the internet.  I have a couple of questions before I proceed:

Would creating a local tile cache from ArcGIS Online Basemaps violate any usage terms.  Without doing any research (yet), I am guessing the local tile cache would be subject to the same usage restrictions of the ArcGIS Online basemaps?
Assuming this is legal, is there anything in place that would prevent me from doing so?  Does anyone know if the ArcGIS Online servers start blocking an IP after n tile requests in n time period?

Any other thoughts,comments,suggestions welcome.


texts: Which URL to enter to access Geoserver via a browser?

texts:
I installed Geoserver using ssh on a server which is hosting a website that I maintain. What will be the URL I need in order to access Geoserver via a browser?

I have tried using the the server host IP. https://nnn.nn.n.nnn/geoserver.

Edit. The server is LAMP and doesn't have Tomcat. I installed the Linux binary and have run the sh. start.sh from the bin directory.


texts: QGIS Methods for Batch/multiple shortest path calculations

texts:
Is there an established method for running multiple shortest path analysis based on multiple start and end points?

I have a model outputting coordinate data and a road network with varying road types. I need to be able to run multiple  shortest path calcs and have those routes plotted.

I am fairly new to qgis and python.


texts: Restoring public.layer_styles table from PostGIS dump

texts:
I created a PostGIS database on a PostgreSQL server and stored some QGIS styles (in public.layer_styles). When dumping and restoring the database only this table gives me an error and stays empty.

I tried only dumping the public schema as .backup and .sql. What's wrong with xml?

 pg_restore: [archiver (db)] COPY failed for table "layer_styles": FEHLER:  invalid XML content
DETAIL:  line 1: StartTag: invalid element name
<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
 ^
CONTEXT:  COPY layer_styles, Zeile 1, Spalte styleqml: â€ž<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="2.8.1-Wien" minimumScale="...â€œ

-- Dumped from database version 9.3.10 to version 9.4.5 <- is this the problem?


texts: Vector Tiles Projection - PostgreSQL/GeoServer/OpenLayers(Mapbox GL)

texts:
I have a dataset in PostgreSQL that is in EPSG:27700, this is read by GoeServer and converted to .pbf files (using GeoWebCache). I have tile gridsets enabled for ESPG:27700,ESPG:4326,ESPG:900913. I am then rendering the tiles using OpenLayers/Mapbox GL.
When using ol.map I get tiles returned and all is working
var map = new ol.Map({
    target: 'map',
    view: new ol.View({
      Projection: 'EPSG:4326',
      center: schladmingWebMercator,
      zoom: 18
    }),
    layers: [new ol.layer.VectorTile({
      source: new ol.source.VectorTile({
        tilePixelRatio: 1, // oversampling when > 1
        tileGrid: ol.tilegrid.createXYZ({maxZoom: 19}),
        format: new ol.format.MVT(),
        defaultDataProjection: 'EPSG:4326',
        url: 'http://someserversomewhere:8080/geoserver/gwc/service/tms/1.0.0/' + layer +
            '@EPSG%3A'+projection_epsg_no+'@pbf/{z}/{x}/{-y}.pbf'
      })
    })]
  });

However changing to mapbox.gl.map it no longer returns tiles as the map centre is outside my extents
var map = new mapboxgl.Map({
        container: 'map',
        minZoom: 15,
        maxZoom: 18,
        style: customStyleJson,
        center: centrepoin,
        zoom: 15,
        url: 'http://someserversomewhere:8080/geoserver/gwc/service/tms/1.0.0/' + layer +
            '@EPSG%3A'+projection_epsg_no+'@pbf/{z}/{x}/{-y}.pbf'
    });

I'm certain the projection is incorrect or not converting as required but I'm struggling to understand where in the connection a change needs to be made. Could somebody explain at what point a different projection can be requested (without any permanent transformation of data)? I notice the projection setting to 4326 in ol.map, is this possible in Mapbox GL?


texts:
I'm trying to get the distance between two points with JTS.
The problem is that this function with JTS

new Coordinate(Lon1,Lat1).distance(new Coordinate(Lon2,Lat2))

Is giving a different result from this function used in postgis(wich provides the good distance)

ST_Distance_Sphere(point1,point2)


texts: Splitting vector file by date (year) but column containing dates in Day/Month/Year format?

texts: I am using national accident data in QGIS and quite frankly the dataset is huge as it contains every recorded accident between 2005 and 2014. I've split datasets by field for a number of vector layers. The issue i'm having is splitting the accident data by year as the date is contained as a field in the DD/MM/YYYY format. I've tried a number of things but as the dataset is so large it typically crashes QGIS.


texts: Cannot load QGIS Processing Plugin due to an error when calling its classFactory() method

texts:
Python error: Couldn't load plugin 'processing' due to an error when calling its classFactory() method See message log (Python Error) for more details.
I've tried redownloading QGIS multiple times, but I am unfamiliar with python so I do not know where the error here is.
I'm using QGIS 3.16.0, but the error still occurred when I downloaded 3.10.0
Couldn't load plugin 'processing' due to an error when calling its classFactory() method
OSError: [WinError 193] %1 is not a valid Win32 application
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\utils.py", line 334, in _startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python/plugins\processing\__init__.py", line 50, in classFactory
    from processing.ProcessingPlugin import ProcessingPlugin
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python/plugins\processing\ProcessingPlugin.py", line 47, in
    from processing.gui.HistoryDialog import HistoryDialog
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python/plugins\processing\gui\HistoryDialog.py", line 33, in
    from processing.gui import TestTools
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python/plugins\processing\gui\TestTools.py", line 34, in
    from numpy import nan_to_num
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\Users\Sophia Weng\AppData\Roaming\Python\Python37\site-packages\numpy\__init__.py", line 140, in
    from . import _distributor_init
  File "C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\Users\Sophia Weng\AppData\Roaming\Python\Python37\site-packages\numpy\_distributor_init.py", line 26, in
    WinDLL(os.path.abspath(filename))
  File "C:\PROGRA~1\QGIS3~1.16\apps\Python37\lib\ctypes\__init__.py", line 356, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: [WinError 193] %1 is not a valid Win32 application


texts: Shapefile created in QGIS will not display properly in ArcGIS

texts:
I was asked to created a point shapefile of a cemetery. The points display correctly in QGIS but not in ArcGIS. I think the issue is with the coordinate system I created it in (WGS 84 / UTM zone 17N).

I've linked the file below for you to try:

https://www.dropbox.com/sh/lv7dc43elfc031r/AADtbioR1wGzmTLikMCKnqOwa?dl=0

What mistake am I making? There is something I am doing wrong during the creation of the shapefile in QGIS.

Edit1: I originally gathered the data using a Garmin GPS unit tagging a sample size of spots. I uploaded the file onto QGIS and created the polyline and polygon shapefiles you can see in the example. These display fine in ArcGIS Desktop.
First, let me be clear that I am new to GIS and even newer to QGIS (I completely self learned the application). I created the point file from scratch by referring to the satellite image plugin of Google and Bing Maps. I also referred to a paper "map" of the cemetery from the 1950's that I won't upload for privacy reasons. My client just wants a digital map to reference and print when needed. They are not GIS professionals but clerks who want to be able to quickly go onto ArcGIS, see if the plot is free or not. Back on topic.

Looking back while I was creating the points, I believe I did something that mkennedy is referring to - changed the geographic coordinate reference. Big rookie mistake on my part.

The thing is, I'm confused and frustrated. Very frustrated. I don't have access to ArcGIS at the moment and cannot play around.


texts: Speeding up raster to point set function

texts:
I have written the following function that takes a raster (always non-skewed) and returns a set of points. As I'm executing it millions of times in my queries, I'd like for it to be super-optimized (my queries currently take hours to complete).

Question: Can the code below be improved?

Who knows a trick for doing this operation faster? Also if the below is not correct, don't be shy to point that out either. Oh, and if you think this is actually a half-decent way of doing it, don't be shy to point that out either :D

For my input (average cells: 38 x 21), the below runs at approx 15ms per raster (rasters are in memory).

Platform is PostGIS 2.0/PostgreSQL 9.2.4.

Thank you.

CREATE OR REPLACE FUNCTION CVL_RasterToPoints
    (
        the_raster raster,
        OUT pt geometry
    ) RETURNS SETOF geometry AS
    $$
    SELECT
        ST_Translate(
            PT.pt,
            (IDX.i - 1) * ST_PixelWidth($1),
            (IDX.j - 1) * ST_PixelHeight($1)
        )
    FROM
        (SELECT
            generate_series(1, ST_Width($1)) AS i,
            generate_series(1, ST_Height($1)) AS j
        ) IDX,
        (SELECT
        ST_SetSrid(
                ST_Point(
                    ST_UpperLeftX($1) + (ST_PixelWidth($1) / 2),
                    ST_UpperLeftY($1) + (ST_PixelHeight($1) / 2)
                ),
                ST_SRID($1)
            ) as pt
        ) PT
    WHERE
        ST_Value($1,IDX.i,IDX.j) = 1
    $$ LANGUAGE sql IMMUTABLE STRICT;


texts:
I'm trying to use EPSG:4479 in QGIS.

1) If I search for EPSG:4479, it returns a "Imported from GDAL" EPSG:4479, but it does not work. At the console, the following message is displayed:

ERROR 6: EPSG PCS/GCS code 4479 not found in EPSG support files.  Is this a valid
EPSG coordinate system?

It should not be a problem, because I can always add a new CRS. But...

2) I tried to define a new CRS in QGIS 2.2, without success.

I go to Settings -> Custom CRS... and fill the new CRS with:
Name: China Geodetic Coordinate System 2000
Parameters: +proj=geocent +ellps=GRS80 +units=m +no_defs

These parameters are copy and pasted from EPSG:4479.

I can use the calculator at the botton to test the coordinate transformation, and it works.

But when I try to save the new CRS, it complains with the message: The proj4 definition is not valid.

What am I doing wrong?

There is a quite similar question around, but not answered. But I'm not using a custom proj4 library.

I'm using QGIS 2.2, PROJ.4 Version 470, on top of the OSGeo Live 7.9.

The same proj4 string works in PostGIS. I am able to do:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 4479, 'EPSG', 4479, '+proj=geocent +ellps=GRS80 +units=m +no_defs ', 'GEOCCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Geocentric X",OTHER],AXIS["Geocentric Y",OTHER],AXIS["Geocentric Z",NORTH],AUTHORITY["EPSG","4479"]]');

select ST_AsText(ST_Transform(ST_SetSRID(ST_Makepoint(114,30), 4326), 4479));

It returns the expected result:

"POINT(-2248544.5475367 5050313.74152924)"


texts: Multiple OL views for single OL Map instance

texts: I need to show multiple views of the same OpenLayer Map instance. Each view would have a different extent, but display the same layers visible on the Map instance. Is that possible? It's important to share the Map state between both views, so I don't believe creating separate Map instances is viable.


texts: Geocortex Essentials: Read data/web map I don't own from a shared private group on AGOL

texts:
I have data and a web map created and hosted in Org B and shared with me via a private group. I am in Org A.

So I was able to successfully add the shared web map as a site reference, but the data (hosted in Org B) was added with no errors or warnings but it does not render on the map in Manager. I opened the viewer to investigate further and I got a message saying it could not load the layers in question.

I tried to add the layers (hosted in Org B, but shared via a group to Org A, and using the URL from Org A) to the map in Manager directly from their service URL, but GE gives me no errors or warnings again in the Manager and just gets hung on the "Add Map Service" window.

Is it possible for GE to consume data hosted from a separate organization and shared via a private group without linking the Manager to multiple AGOL Organizations?


texts: How to use v.clean in QGIS to fix polygon topology

texts:
I have QGIS 3.18 Zurich with GRASS 7.8.5 installed on my 64bit Windows machine. I have digitized 1600+ polygon boundaries that now have sliver polygons, gaps, and overlaps which I've confirmed via the Topology Checker Plugin.
Data are in EPSG 32629 (WGS84 UTM 29N). I've created a key variable FID which was calc'd as @row_number. I am attempting to run v.clean on this data in QGIS with the following parameters:
Input Feature Type = area
Cleaning tools = bpol, rmarea, rmline, rmsa, rmdupl, rmdac
Thresholds = 5,5,5,5,5,5  (this is 5m; I have also tried no thresholds)
Output type = auto
Combine tools with recommended follow-up tools = yes
V.in.ogr snap tolerance = 0.05 (5cm)
v.in.ogr min area = 10.0 (10 sqm, smallest polygon is 11sqm)

However I keep getting blank "Error" and "Cleaned" out puts and receive the following errors in the log:
WARNING: Number of centroids exceeds number of areas: 1619 > 1448
WARNING: Number of incorrect boundaries: 202
WARNING: Number of centroids outside area: 171
ERROR: Key column <> not found
WARNING: Unable to determine input map's vector feature type(s).
WARNING: No attribute table found -> using only category numbers as attributes
WARNING: Output layer is empty, no features written

What am I doing wrong?
Data are here: https://drive.google.com/file/d/1lYIVkftS_XLRF1CQcNsjNI_s5ToegC8w/view?usp=sharing


texts:
I started off with a SpatialPointsDataFrame, containing locations from animals wearing GPS loggers. Aside from the coordinates, the data frame also contains the animals' ID and a timestamp.

From this SpatialPointsDataFrame, I created lines that show the trajectory of each animal over the tracking period. I've used the following code:

spdf4tri#the SpatialPointsDataFrame

ids=unique(spdf4tri$id)
trajectory <- list()
 for (i in ids){
   spdf4bird<-subset(spdf4tri, id==i)
   birdtrajectory<-SpatialLines(list(Lines(list(Line(spdf4bird)), "id")))
   trajectory[[i]]<-birdtrajectory
   print(i)
 }

Now, I have a problem, because trajectory is of the class 'list'. I cannot figure out how to unlist it.


texts:
I have some line layers ,I wan't to create 0.5km and 1km buffer around the lines.

How to create the the buffer in OpenLayers 3 and which layer use like wms,wps or wfs?


texts: How to calculate a percentage of an area covered by polygons ArcGIS 10.1

texts: I created an Aspect Layer with 9 Zones on which I want to use the Tabulate Area Tool with another polygon layer namely Woods. So what I actually need to know is the percentage, covered by woods for each of these areas. I tried it with tabulate area but without success, the intersect tool might be useful here but I can choose the polygon layer there.


texts: Can QGIS export Geospatial PDF for Avenza PDF map app?

texts: Since QGIS uses GDAL, does any map exported as a PDF become a Geospatial PDF? I am hoping to import maps from QGIS into the Avenza PDF Map app on my phone. Previous posts mentioned exporting as Geospatial PDF's and GeoTiff's, but none had an exact answer or process.


texts: Geoserver WMS distortion when reprojecting on high zoom

texts:
I have a GeoServer WMS that I load in a Leaflet map (keeping the map default CRS 3857) as a L.tileLayer.wms. The problem is that when I zoom in on this tile layer, the pixels get deformed and appear split, as shown below in the image (First img: Leaflet-map CRS: 3857; Second img: QGIS project CRS 3857). The second image is my desire result.

(Leaflet)

(QGIS)
Minimal code needed for a reproducible example:
var map = new L.Map('leaflet', {
    center: [28.3300,-16.490740],
    zoom: 15
});

L.tileLayer.wms('http://sigred.oapn.es/geoserver/SistemaREMOTE_modis/wms?',{
    layers: "SistemaREMOTE_modis:ref_prod_media_resumen",
    transparent: true,
    format: 'image/png'
}).addTo(map);

I tried changing the CRS of the Leaflet map to 4326, but that didn't work, and also that way my OSM basemap won't place right.
Also If I visualize the WMS with OpenLayers from GeoServer it looks correctly (http://sigred.oapn.es/geoserver/SistemaREMOTE_modis/wms?service=WMS&version=1.1.0&request=GetMap&layers=SistemaREMOTE_modis:ref_prod_media_resumen&styles=&bbox=-965199.01397601,3190740.34917998,834315.241888609,4798237.84077302&width=768&height=686&srs=EPSG:25830&format=application/openlayers)
What I am doing wrong? Is there a way to solve this problem?


texts: Appending a table asset with ImageCollection pixel values for specific dates and coordinates in Earth Engine

texts:
In Earth Engine I have a table asset of coordinates and dates (properties `LAT', 'LONG', 'DATE') for which I would like to return a pixel value from my ImageCollection dswe_ic, which gives a monthly DSWE (Dynamic Surface Water Extent) value.

// get a DSWE value (integers of 0, 1, 2, 3, 4, or 9)
var dswe_ic = ee.ImageCollection.fromImages(months.map(function(m){
  var start = start_date.advance(m, 'month'); // Advance start date by m months
  var end = start.advance(1, 'month');        // End date is start month +1
  return img_indices_all.filterDate(start, end)
          .select(['dswe'])
          .reduce(ee.Reducer.min()) //  Produces band 'dswe_min'
          .set('date', start.format("YYYYMM")) // Sets yr/mo as property 'date'
          .cast({'dswe': 'uint8'})  // Adds empty band 'dswe'
          .remap([0,1,2,3,4,8,9], [0,1,2,3,4,9,9]).rename('dswe'); // remap takes the 1st band unless otherwise specified
          // return a monthly DSWE value, taking the lowest DSWE value (i.e. highest confidence value)
          // of the month if there are multiple DSWE values to pick from
          // (due to potential availability of more than 1 high quality Landsat scenes)

I can map over my coordinates of interest with dswe_ic easily enough and generate a CSV list of monthly pixel values from my start date (May 21, 2003) to my end date (Nov 17, 2016).

var results = points_of_interest.map(function(f) {
  return dswe_ic.map(function(i) {
    var mean = i.reduceRegion({
      geometry: f.geometry(),
      scale: 30,
      crs: 'EPSG:4326',
      reducer: ee.Reducer.first(),
    });

    return f.setMulti(mean).set('date', i.get('date'))
  })
})
results = results.flatten()

Export.table.toDrive({
  collection: results,
  description: 'results_all_coords_allmonths_from_may2003_to_nov2016,
  fileFormat: 'CSV'
});

However, I am more interested in directly appending my original table with a new column for the pixel value from dswe_ic for each particular 'DATE' and coordinate.  Unfortunately, my GEE skills are quite limited to do this.  Any suggestions?


texts: Borders appear when reducing the global opacity of multiple raster images in QGIS 3.2.3

texts: I want to reduce the opacity of multiple raster images and when I use the global opacity slider to do so, black borders appear between my images. I would want to eliminate these borders.


texts: Writing plugin to develop composer attribute tables?

texts:
I have no experience of plugin development, but I'm wondering if it would be possible to write a plugin that solved the issue of multi-line tables in print composer http://hub.qgis.org/issues/10273

I now have several GIS projects in QGIS, having found this software only a few months ago, after having several teething troubles answered in this site, the only thing preventing me from just printing up my projects is the lack of multi-line support in print composer tables. I'm therefore looking at every possibility for solving this problem, including programming it myself, but I don't want to start the process if it's not even possible.

Does anyone know if such a plugin is possible, or that just not the sort of thing that plugins can do?

It seems such a small thing.


texts: Overriding method of qgis object does not work

texts:
I try to reimplement the renderDiagram() method from QgsHistogramDiagram to add some labels directly to the diagram (using the position argument).
So far, I tried to set up my own class to override the method:

class OwnHistogram(QgsHistogramDiagram):
    def __init__(self):
        super(OwnHistogram, self).__init__()

    def renderDiagram(self, feature, renderContext, diagramSettings, position):
        # Here I want to see if function is entered:
        print('Try to override')
        super(OwnHistogram, self).renderDiagram(
            feature, renderContext, diagramSettings, position)

But the function is not entered, while Histogram is shown as if nothing has changed.
Is it not possible to get into the C++-wrapped methods, or did I miss something?


texts: Adding Microsoft Access Table to ArcMap?

texts:
I am looking to map gps coordinates from a Microsoft Access table (.mdb) in ArcMap. I successfully connected ArcMap to Microsoft Access through an OLE DB connection, and can see the Access table in both ArcCatalog and the Add Data dialog box. However, when I attempt to add the table I receive the following warning:

  A selected item could not be added to the map. General Function
  Failure.

After searching online for a solution I thought maybe ArcMap could not read the file because of an incompatible driver. I am running Microsoft office 2016 (64-bit) on a 64-bit machine. However I read that since ArcGIS is a 32-bit program the correct driver for this setup would be Microsoft Office Access Database Engine 2007, which is the driver currently installed on my computer.

All other documentation I have found online indicates that this should be a fairly simple process, and I can not find any mention of the problem I am having.

Any suggestions on why this is happening?


texts:
I am working on a case now, in which the user can select 2 or more polygons on the map. I want to merge these polygons each time user clicks on a new one. This is what I did so far ... But seems it doesn't work well because I need to clear the previous layer before calling it again .

$("#merge").click(function() {
   console.log('Merge Request !');
   var merged;
   var selected = [];
   geojson.on('click', function(e) {
      selected.push(e.layer.toGeoJSON());
      if (selected.length == 2){
         var pol = turf.union(selected[0], selected[1]);
         var merged = L.geoJson(pol,{'style': function() {
               return {
                  'color': 'white',
               }
            }
         }).addTo(map);
         merged.bringToFront();

         selected = [];
         selected.push(pol);
      }
   });
});


texts: How to restore default values in Mapserver?

texts:
I'm using Mapscript reading a configuration of my map file, so i set some values in my mapObj. But in certain point of my code i need to restore the default values locate in map file. I tried to use the method "applyConfigOptions()", but didn't work.
Is there a way to do that?


texts:
When I use the merge tool on this polygon mesh, I end up with a multipart polygon with large gaps.

Before this, I was getting an error when trying to run a merge or dissolve because of "Polygons with overlapping Z values", so I exported a copy of the feature class without Z values to run the merge on.


texts: Can we set the latlng start point : leaflet

texts:
I am working with the leaflet api.

I am trying to figure out how to change or set the latlng start point in leaflet map.

  Currently i am working on use image as map,,,,the latlng(0,0) is at
  the top left of the image or container of that image...

  If its possible that i change it to left-bottom

Snapshot of map :

If its possible to change the position of origin (0,0) from top-left to bottom-left,,,,


texts: Shapefile to GeoServer publishing to store error

texts:
I am trying to publish a shapefile (archived) to GeoServer, by using a Python script, then display the shapefile via WMS in JOSM.
When I publish it to a new Workspace it works fine, it is displayed in JOSM. When I publish it to an existing Store within the workspace or when I create a new store within the workspace, I cannot retrieve it via WMS in JOSM. It is only displaying the first shapefile sent to that workspace.

Code is below:

from geoserver.catalog import Catalog
from geoserver.util import shapefile_and_friends

geoserver_uri = 'http://localhost:8080/geoserver/'
from geoserver.workspace import ResourceInfo
cat = Catalog(geoserver_uri + 'rest')
# topp = cat.get_workspace("topp")
# shapefile_plus_sidecars = shapefile_and_friends("Tasmania roads")
# print(shapefile_plus_sidecars)
wks_name = 'newWorkspace'
uri = geoserver_uri + wks_name
store_name = 'cityName'
data = r"C:\Users\this_pc\Documents\cities\cityName.zip"

def check_workspace(workspace_name):
    check_workspace = False #create a variable to test if workspace exists already
    if(cat.get_workspaces(workspace_name)):
        check_workspace = True #set the variable to true, if workspace exists
    return check_workspace

#check_workspace(wks_name)
#check_store = False

def create_data(store_name, data, data_name, wks_name):
    cat.create_featurestore(store_name, data, workspace=wks_name)
    cat.add_data_to_store(store=store_name, name=data_name, data=data, workspace=wks_name)

if check_workspace(wks_name): # if workspace exists, create the featurestore
    print("Workspace exists")
    if(cat.get_store(store_name, workspace=wks_name)):
        cat.add_data_to_store(store=store_name, name=store_name, data=data, workspace=wks_name)
        print("Store found. Data has been written.")
    else:
        print("Store does not exist. Store has been created. Data has been written to the store.")
        create_data(store_name, data, "Test name", wks_name)

else: #if workspace does not exists, create the workspace and the featurestore
    cat.create_workspace(wks_name, uri)
    create_data(store_name, data, "Test name", wks_name)
    print("Workspace does not exist. Workspace and store has been created. Data has been written to the store.")

Edit: Also, when I publish it, the zip file is erased from local drive.

  Traceback (most recent call last):
    File "C:/Users/paul.nicoara/Documents/Slack/Mapillary_DB/save_geoserver.py", line 40, in
      create_data(store_name, data, "Test name", wks_name)
    File "C:/Users/paul.nicoara/Documents/Slack/Mapillary_DB/save_geoserver.py", line 27, in create_data
      cat.add_data_to_store(store=store_name, name=data_name, data=data, workspace=wks_name)
    File "C:\Users\paul.nicoara\PycharmProjects\Mapillary_DB\venv\lib\site-packages\geoserver\catalog.py", line 429, in add_data_to_store
      with open(bundle, "rb") as f:
  FileNotFoundError: [Errno 2] No such file or directory: 'C:\Users\this_pc\Documents\cities\Bandung.zip'


texts: Clustering geographical data based on point location and associated point values

texts:
Given data points with longitude, latitude, and a third property value of this point. How can I cluster points into groups (geographical sub-regions) based on the property value? I searched by google and figured out that this problem seems to be called "spatial constrained clustering" or "regionalizing". However, I am not familiar with handling geographical data and haven't get an idea about what kind of algorithms are good, and which python/R packages are good for this task.

To give a more intuitive idea about what I want, let's say my data scatter plots are as following:

So each dot is a point, x is longitude, y is latitude, and colormap shows whether the value is big or small. I want to divide those points into sub regions/groups/clusters based on location and similarity of values. Like the following (it is not exactly what I want, just to show a intuitive idea.):

So how can I achieve this?


texts:
I have three tables in Postgis:

points called location(gid (pkey),name) example (1,Paris)
temperatures named Temp(idTemp(pkey),TempValue,date) example (1,12,2015-01-01)
relation between the two: TemperatureAtLocation(gid,idTemp)

The aim is to display all temperatures detected in Paris in Qgis if possible when hovering or clicking on the city symbol or any user friendly solution.
I checked if any actions available for vector layer exist but got little success.


texts: Visualisation of spatio-temporal data set?

texts:
I'm investigating visualisation of a spatio-temporal data-set.

What examples are available, online or offline, of high-quality spatio-temporal data visualisation?


texts: Calculating slope for points with downhill slopes rendered as negatives in ArcGIS

texts:
I have a DEM and a series of points. I'd like to measure the change in slope along the series of points. Distinguishing between uphill and downhill slopes matter to me, and ArcGIS' slope function renders all slopes as absolute values.

I'm working with ArcGIS' Spatial Analyst and 3D Analyst toolkits.

This is what I've done so far.

Created points
Run the "interpolate shape" tool
Run the "slope" tool
Run the "extract multi values to points" tool to add two fields Elevation and Slope to the points feature class.

I have followed the advice here (http://blogs.esri.com/esri/supportcenter/2012/03/09/calculating-the-slope-of-a-hiking-trail/) to change the downhill slope values into negatives, but have not succeeded. When I run the code (below), nothing happens.

fc = r"C:\Points\Test.gdb\Points_3d"
cursor = arcpy.da.UpdateCursor (fc, ["Elevation","Slope"])

firstRun = True

for row in cursor:
   if firstRun:
      oldValue = row[0]
      firstRun = False
   else:
      if oldValue < row[0]:
         row[1] = row[1]
      else:
         row[1] = (-1)*row[1]
      oldValue = row[0]
      cursor.updateRow(row)
del row
del cursor


texts: Selecting feature and highlighting based on feature properties in OpenLayers 3?

texts:
I need to select link using link properties, I am able to find the feature but not able to highlight programmatically. please help..

In openlayers 2 we have selectFeature but not sure what is in ol3


texts: QGIS Print composer's Atlas export file options missing Tiff on Mac

texts:
Trying to create a workflow using QGIS for tree canopy inventory projects that works in both Windows and Mac.
I do not have a Mac, but doing peer review is on Mac.

There seems to be an issue with the Print Composer's Atlas export image not having the option to select tiff for the export output.
Screen shots are here and below.

Can anyone here with Mac verify the lack of export to tiff option or post screenshot of where the option is located?


texts: In mapbasic,how to register a table in xlsx,not in xls?

texts:
I have tow files that only differ in  extention name, one xls ,the other xlsx.both contained the same things,same structure.there's no problem when I registered a table in xls,but when in xlsx,an error occured.

code list:

Include "mapbasic.def"
Include "menu.def"
Include "icons.def"

Declare Sub main
Declare Sub Combine

Sub main
   Register Table "C:\Documents and Settings\Administrator\desktop\test\fdd.xlsx"
   TYPE xls

End Sub


texts: Converting high precision dataset in file geodatabase to low precision using ArcObjects?

texts:
A high precision dataset was loaded from a file geodatabase into arcsde 9.3.1. The dataload was achieved by copy/paste in ArcCatalog.

A custom tool (arcobjects) that extracts from the dataset to a low precision personal geodatabase is now failing when attempting to create the output feature class.

It is not feasible to upgrade the personal geodatabase to high precision since all the other datasets it extracts from sde are low precision.

Can the high precision dataset be downgraded?


texts: How to compute azimuth between QgsPointV2s?

texts:
The QgsPoint class has an azimuth() method but the new QgsPointV2 class does not. What's the best way to compute the azimuth between QgsPointV2s? I couldn't find a convenient way to convert V2 to the old version, so my fallback is to construct a QgsPoint from the V2's x and y:

self.traj = QgsLineStringV2(...)
p_i = self.traj.pointN(self.i) # this returns a QgsPointV2
p_j = self.traj.pointN(self.j)
a = QgsPoint(p_i.x(),p_i.y()).azimuth(QgsPoint(p_j.x(),p_j.y()))

Is there a better way?


texts: Using ST_Union to merge tables containing lines?

texts:
I have to build a street network with PostGIS that will be used with pgRouting.

I have two tables:

the_network (182 lines): is a MULTILINESTRING table. Having columns "road_name" (varchar),  "road_length" (numeric)
buildings (104 lines): is a POINT table. Having "building_id" (integer)

The following code creates lines between buildings and the network. The 104 created lines are stored in "buildings_to_network_lines_table" which is GEOMETRY type.

SELECT ST_ShortestLine(a.geom, ST_Collect(b.geom)) AS geom,
       ST_Distance(a.geom, ST_Collect(b.geom)) AS shape_length INTO buildings_to_network_lines_table
FROM  buildings a, the_network b
WHERE ST_Dwithin(a.geom, b.geom, 1000)
GROUP BY by a.gid, a.geom;

Then I just want to create one table named the_network_full that contains all the lines of the_network and buildings_to_network_lines_table tables. I've tried to achieve this by a ST_UNION() function:

SELECT ST_Union(a.geom, b.geom) AS the_network_full
FROM  the_network a, buildings_to_network_lines_table b;

The results look good but unfortunately the resulting "the_network_full" table (type GEOMETRY) contains 18'928 lines (182*104) instead of 182+104=286 lines.

How can I merge these two tables into a new one?


texts:
I have a road network shapefile.

After creating Network dataset I get a graph with vertices and edges.
I need to create a tool for ArcGIS (through a Python code) which enables to calculate depth of all vertices in a network according to the formula:

Number of neighbours one step away multiplied by number of steps + neighbours two step away multiplied by number of steps + three steps away...

For example see image: red point on the picture has 3 neighbours one step away, 2 neighbour 2 steps away. So the depth for the red point will be: 3x1+2x2=7.

I need to calculate this for all vertices in my network. It's quite easy but unfortunately I don't have any experience with Python.


texts: Improve display of fine resolution raster

texts:
In the above map, I am showing yields on some parts in Iowa at a pixel resolution of 100m. However, in this case due to the relatively small number of pixels (and high resolution), the map is not very clear. Is there any way I can smoothen out the raster so that it displays nicely using ArcGIS?

Thanks to @ Michael Miles-Stimson for suggestions, the yields that I am displaying are the 'VALUE' of the raster. Therefore it is not possible to do something like a focal mean (which destroys the raster). I would rather avoid changing the attribute table to use the yields as 'VALUE' because there are several other variables I want to plot as well, and this would loose information and be very time-consuming for my raster.


texts:
I'm working with 19th century land tax records that provide a direction (expressed as NW, SE, etc.) and a distance in miles from the court house. The info in the Land Tax Records (LTR) is in effect a version of polar radiation. For each reading you start at the court house as your 0,0. The records have: LastName; Direction, Distance or LastName, Distance, Direction. Such as Smith, 4, SW that translates to Smith, 4 miles southwest of the court house.

Is there a way to place a symbol/marker/whatever at that point using QGIS 2.0? The end result should roughly result in the shape of the county in 1860.

It is possible and very cumbersome to do that with CAD by using the direction and distance alignment on a line and then placing a loci at the end of it. I can use the direction and distance measures by placing them into Excel, plotting them in Vectorworks, and then exporting them to QGIS. I'm trying to cut out the intermediate steps by having QGIS do the work. I can convert the directions to azimuths in Excel and the distances to whatever I need. But having 400 separate lines for a county is cumbersome if there's a more elegant solution.


texts: Changing a markers location at set intervals

texts:
I am trying to plot coordinates of a marker using a GeoJSON feature. I am able to plot my points fine but my issue is, when the interval happens the old marker stays on my leaflet map rather than removing the old marker or changing the coordinates of the original marker.

The code below is in my js file which has a set interval to read from a URL with my coordinates and this is where my marker is added to my map.

    window.setInterval(function nameSat()
    {
 $.getJSON('http://10.175.129.36/meo_coordinates_v2.json', function(data) {
         lat =`${data[1].Latitude  }`
         lon =`${data[1].Longitude }`

         //console.log("Lat: "+ lat + "Long: " + lon);

 geojsonFeature = {
    "type": "Feature",
    "properties": {
        "name": "Satellite",
        "marker-color": "#ff0080",
        "marker-size": "small",
         "popupContent": "O3B M003
Longitude: "+ lat + "
Latitude: "+ lon
    },
    "geometry": {
        "type": "Point",
        "coordinates": [lon, lat]
    }
};

console.log("Lat: "+ lat + "Long: " + lon);

    var meo =L.geoJSON(geojsonFeature , {
    onEachFeature: onEachFeature
});
meo.addTo(mymap);

 })

    }, 10000);


texts: Defining polygon object by identifying lines in R

texts:
I have a dataset contains lines and I have imported them into R. I want to take a close look at the coordinates of them and define the identical first and last coordinate of each point if there is some -->(looking for polygon). Therefore, I am using Slot which makes me able to have a close look at the details of the desired object. My final goal is to define the number of identical point coordinates(First and last) for each line in order to discover the number of the potential polygon in my data.

How many of lines objects have potential to be a polygon?

To do so, I have done several steps: In the first step, I read my data into R. The second step, I have used slot to have a close look at the coordinate of each point(Sequence of points represent line object). the third step: I have tried to define the number of identical points but I have faced with an error says CRDs not found

at the following, you can tack a look at the codes

enter library(maptools)
#Read data directly from National Geophysical Data Center (NGDC) coastline
#extractor.
shorelinedat="http://www.asdar-book.org/RC1/datasets/auckland_mapgen.dat"
#Assign CRS
llCRS <- CRS("+proj=longlat +ellps=WGS84")
#Read data from mapgen into a SpatialLines object.
auck_shore <- MapGen2SL("auckland_mapgen.dat", llCRS)
#Required code to identify the lines.
lns <- slot(auck_shore, "lines")
table(sapply(lns, function(x) length(slot(x, "Lines"))))

Here is the code in which I faced with the error

#identifying the number of identical coordinates
islands_auck <- sapply(lns, function(x) {
+ crds <- slot(slot(x, "Lines")[[1]], "coords")
+ identical(crds[1, ], crds[nrow(crds), ])
+ })

This is the error

Error in +crds <- slot(slot(x, "Lines")[[1]], "coords") :
object 'crds' not found


texts: Overlaying Landsat 8 (raster) and shapefile using ArcMap

texts:
I would like to know what I can do correct such problem, those images are not overlaying while they have the same projection system.
I used even the extent of the shapefile to change the one for the Landsat but the problem is still available.


texts: How to run special characters in field calculator?

texts:
How to run a special character in field calculator in arcmap. Special character is not working in arcmap field calculator.

i need Like this (Três into Tres)

I am using this code in python

!ADDRESS_N!.replace("Três","Tres")