texts: Looping workflow across attribute table using Google Earth Engine

texts:
I have a shapefile asset loaded with an attribute column called "OBJECTID". I would like to run the following script for every ID number in the attribute table, and produce an output table.

/// define a geometry called Geometry.

//visualize the area you are working with
Map.addLayer(geometry)

// Landat 5 surface reflection data
var L5coll = ee.ImageCollection("LANDSAT/LT05/C01/T1_TOA")
//Filter the cloud cover
.filter(ee.Filter.lt('CLOUD_COVER',40))
//select the desired NIR and R bands
.select(['B3', 'B4'])
//use only the geometry specified
.filterBounds(geometry)

// Landat 7 surface reflection data, fill in the gaps! See USGS pages for more info
var L7coll = ee.ImageCollection("LANDSAT/LE07/C01/T1_TOA")
//filter the cloud cover
.filter(ee.Filter.lt('CLOUD_COVER',40))
//select the NIR and R bands
.select(['B3', 'B4'])
//Use geometry specified
.filterBounds(geometry)
//creaate a variable that takes the focal mean of the pixels surrounding the missing gaps of L7 imagery
//the variable is then used to return a blended image which contains filled in gaps using .blend func.
.map(function(image){
  var filled1a = image.focal_mean(2, 'square', 'pixels', 1)
  // fills in gaps via blending
  return filled1a.blend(image);
})

// Landat 8 surface reflection data, rename the band names. See USGS pages for more info
var L8coll = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
//filter for cloud cover
.filter(ee.Filter.lt('CLOUD_COVER',40))
//use desired geometry
.filterBounds(geometry)
//applies the function over the image collection
.map(function(image){
  //renames the bands to match previous landsat imagery
  return image.rename(['B0', 'B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B9', 'B10', 'B11']);
})
//select NIR and R bands, now renamed to match L5 and L7 imagery
.select(['B3', 'B4']);

// merge L5, L7 & L8
var collection_merge = ee.ImageCollection(L5coll.merge(L7coll.merge(L8coll)));
print (collection_merge)
// add NDVI
var addNDVI = function(image) {
  //GEE has a normalized difference function built in
  var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  return image.addBands(ndvi);
};
// create image collection from merged landsat imagery, filter boundaries and dates, select bands
var collection1 = collection_merge
//filters based on AOI
    .filterBounds(geometry)
    //Filters the input dates
    .filterDate('2000-01-01', '2017-12-31')
    //selects NIR and R bands
    .select(['B3', 'B4'])
    //applies NDVI to image collection merge
    .map(addNDVI).select('NDVI');

//Calculating annual NDVI
    var year = ee.List.sequence(2009,2017);
    var year_func = function(y){
      var range = ee.Filter.calendarRange (y, y, 'year');
      return collection1.select('NDVI').filter(range).mean().set ('Year', y)
    };
    var yearwise_ndvi = ee.ImageCollection(year.map(year_func));
    print (yearwise_ndvi);
    Map.addLayer (yearwise_ndvi)

//Creating time-series chart:
var chart = ui.Chart.image.seriesByRegion ({
  imageCollection: yearwise_ndvi,
  regions: geometry,
  reducer: ee.Reducer.mean(),
  scale: 30,
  xProperty: 'Year'
}).setChartType('ScatterChart')
  .setOptions ({title: "Median NDVI 2009 - 2017",
  hAxis: {title: 'Year', format: 'year'}
});

// chart the NDVI time series
//chart by region uses input geometry and yearwise_ndvi
//ee.Reducer.median is taking the median NDVI for each tile charted

print(chart)
//export the time series as a table via console
//in the console the table can be exported from the chart
var projectyear = 2000
var startyear = projectyear - 1
var endyear = projectyear + 1

// Run the loop
for (var yr= projectyear; yr <= endyear; yr = yr + 1) {

  print('This is year ' + yr)
  print(yr)
}


texts:
I am trying to install the QGIS Desktop application on Windows 10.
Upon a successful install, only the OSGeo4W Shell and the Grass GIS 7.8.5 are there on the apps list (within the QGIS folder), the actual QGIS Desktop app is not there.
I get the .exe file from the regular qggis dot come downloads website ("QGIS-OSGeo4W-3.16.5-1-Setup-x86_64").
I've tried with both the 64-bit standalone installers for 3.16 Hanover and (latest) 3.18 Zurich version. This process has worked just fine in the past.
My computer's virus protection program Sophos Endpoint Agent keeps detecting and clearing a file named nircmd.exe from the QGIS/3.1x/bin folder, wondering if that has anything to do with it (says "Threats cleaned up. NirCmd detected at file path..."). I don't know what I'm doing wrong now or what is missing. Please help.


texts: Standalone app does not show the content

texts:
I am trying to show the shapefile in canvas using standalone application. However it does not show the content. Here is my code and I am using PyQt5, Qgis 2.99, PyCharm 3.6

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt5.QtCore import *
from PyQt5.QtGui import *

app = QgsApplication([], True)
path = "C:/Users/DELL/AppData/Roaming/QGIS"
app.setPrefixPath(path, True)
app.initQgis()
canvas = QgsMapCanvas()
title = "PyQGIS Standalone Application Example"
canvas.setWindowTitle(title)
canvas.setCanvasColor(Qt.white)
layer_info = "C:/Users/DELL/Desktop/GIS_CensusTract_poly.shp"
layer = QgsVectorLayer(layer_info, 'MyShp', "ogr")
layer.updateExtents()
QgsProject.instance().addMapLayer(layer)
canvas.setExtent(layer.extent())
canvas.setLayers([layer])
canvas.zoomToFullExtent()
canvas.show()
exitcode = app.exec_()
QgsApplication.exitQgis()
sys.exit(exitcode)


texts: Extracting data.frame from simple features object in R

texts:
Is there an sf-native (i.e. "correct") way of extracting everything except the geometry column from a simple features object? This works

df <- dplyr::select(as.data.frame(sf), -geometry)

but the select( , -geometry) step feels unnecessary. Also, it doesn't remove the geometry attributes.


texts: How to reproduce geodesic transformation in Google Maps API

texts:
I need to show points on country map of Russia, but not in Mercator projection, it should be in usual view, be more real (like we see it everywhere every day). Trying to figure out different projections, transformations and etc, I've created a test page with Google Maps and deogesic polygon with geojson points of the country, move it to equator and got usual form of Russia, that suits me. Google Maps performs all transformations for me, I just need to mark this poligon as geodesic. But problem is there is no API to convert coortinates in the same way and I don't know how to recalculate coords of particular point, which I pick on Mercator map, to get that point be in right place on transformed map.

After some googling I found project http://bramus.github.io/mercator-puzzle-redux/. Author created a function to move poligon with geodesic transformation. The algorythm looks pretty simple: find center of polygon, evaluate heading and distance of each point of paths, get new coordinates according to new center using spherical geometry functions. But calculations isn't precise (and author admits it): when I move country to point (0,0) and return it back, borders doesn't match (screenshots). Even center point of origin country doesn't match with same geographical point on transformed map. Here is an example - https://jsfiddle.net/unnamed777/rg7tp4zL/ . When I drag Russia, Maps API performs all transitions. Green point is center of geojson poligon (calculated by the same algorhitm as Bramus' function moveTo() via center point, heading and distance). And when I move country, new center of that poligon doesn't match with real coordinates of surface.

So, my question is what algorythm do I need to use to transform points from classic Mercator projection to that with minimal distorsion as if Russia is placed on equator? Actually, I want to end up with svg/png image of country with highlighted regions and function, that receives point picked on google maps and translates it relatively to pixels of my image (https://i.imgur.com/3xqB3Sa.png). But for the beginning I want to understand how to translate coordinates as google does with poligons.


texts: The Calculate Field tool causes an error when used in a geoprocessing service

texts:
I wrote a tool that creates a table whith errors and fills an error field.
Then I published it on the server as a geoprocessing service. But when I run it causes an error

But when I comment arcpy.CalculateField_management(temp_table_error, "Error_text", "'Раскорчевка должна быть заполнена'", "PYTHON_9.3") it works great. I don't understand what should I do.
arcpy.env.workspace = "C:\inetpub.........\GosKontrol_temp.gdb"
arcpy.env.overwriteOutput = True
arcpy.CreateTable_management(arcpy.env.workspace, "Excel_Error_"+str(varOrgname))
temp_table_error=arcpy.env.workspace+"\\Excel_Error_"+str(varOrgname)
arcpy.AddField_management(temp_table_error, "ObjectID_vovl", "LONG", 9)
arcpy.AddField_management(temp_table_error, "Error_text", "TEXT", 255)
UchGoskontr_temp="C:\inetpub\......\GosKontrol_temp.gdb\UchGoskontr_temp"
##########################################################################
expr1=""""Kod_vz"=1 AND "Korch_vz" IS NULL AND "Orgname" = """ + str(varOrgname)

arcpy.MakeTableView_management (UchGoskontr_temp, "stateslyr1", expr1)

#
fieldmappings = arcpy.FieldMappings()
fieldmappings.addTable(temp_table_error)
try:
    fieldmap = fieldmappings.getFieldMap(fieldmappings.findFieldMapIndex("ObjectID_vovl"))
    fieldmap.addInputField("stateslyr1", "OBJECTID")
    fieldmappings.replaceFieldMap(fieldmappings.findFieldMapIndex("ObjectID_vovl"), fieldmap)
except:
    pass

arcpy.Append_management("stateslyr1",temp_table_error,"NO_TEST",fieldmappings)

arcpy.CalculateField_management(temp_table_error, "Error_text", "'Раскорчевка должна быть заполнена'", "PYTHON_9.3")

arcpy.Delete_management("stateslyr1")


texts: Adding boolean yes/no field in ArcGIS Desktop?

texts:
How does one add a field to an ArcGIS feature class with a boolean data type? That is an attribute where the allowed value is only one of a pair such as 1/0, Yes/No, On/Off, Present/Not-present, etc.

A boolean data type is not available in the Data Type pick-list.


texts: PostGIS: Perform a Spatial Join only joining the intersecting polygon with the maximum value

texts:
I'd like to perform a spatial join using PostGIS. I've got table_A containing 80.000 entries and table_B containing roughly 1.500.000 entries. The entries in table_B have a column "number".

My goal is to only join the entries from table_B with the maximum value which intersect polygons from table_A.

For example: In the following image I'd like to only join "polygon 3" from table_B to the blue polygon of table_A. Polygon 1 and Polygon 2 should not be joined because they have lower numbers.

I could also use ArcMap (Standard license) but it crashes due to the size of the shapefile (~ 1.8GB).


texts: Setting OpenLayers line labels by angle?

texts:
I want to write labels near linestring features on my map. I calculated but it works wrong.

My layer style function is like following:

layerFeatures = new ol.layer.Vector({
    source: sourceFeatures,
    style: function(feature, resolution){

        let first = feature.getGeometry().getCoordinates()[0];
        let last = feature.getGeometry().getCoordinates()[1];

        let rotation = -Math.atan2((last[1] - first[1]),(last[0] - first[0]));

                    var style = new ol.style.Style({
             fill: new ol.style.Fill({
             color: 'rgba(255,255,255,0.4)'
           }),
           stroke: new ol.style.Stroke({
             color: '#3399CC',
             width: 1.25
           }),
           text: new ol.style.Text({
             font: '14px Verdana',
             fill: new ol.style.Fill({
              color: '#000000'
             }),
             stroke: new ol.style.Stroke({
              color: '#f2f',
              width: 1
            })
           })
        })

        let label = feature.get("properties").name;

        style.getText().setRotation(rotation)
                    style.getText().setText(label)

        return [style];

    }
})

The directions of labels are wrong. It shold be readable direction, but some of them are reversed.

and my working code is here


texts: Making 95% Volume Contour on ArcMap without Hawth's Tools or GME?

texts:
I have some points that represent animal locations over the course of a few months. I am emulating the work of previous researchers in this field. I'm still new to GIS, so I previously thought that I wanted to run my animal locations through a Kernel Density Estimate. However, re-reading the previous literature and comparing the maps, I think I need a particular kind of KDE. The previous research was published in 2011, so the author used Hawth's tools to conduct a Fixed KDE with volume contours (referenced here and below). Unfortunately, Hawth's doesn't work with Arc10 and it's analogue, GME doesn't work with ArcMap 10.6, which is on all the university computers.

Is there another tool that I can download to make this percent volume contour?

Here is a picture of the KDE from the cited report

Management of pinyon-juniper woodlands at Kirtland Air Force Base: Pinyon Jay summer and winter home ranges and habitat use 2009 final report:

  For the kernel analysis we used the Hawth’s Tool, Fixed Kernel Density
  Estimator. Kernel analysis is a nonparametric statistical method for
  estimating probability densities from a set of points. When used to
  analyze home range data, kernel density methods describe the
  probability of finding an animal in any one place. The Fixed Kernel
  Density Estimator calculates a fixed kernel density estimate and
  produces contour lines representing the boundary of the area that
  contains a specified percent of the volume of a probability density
  distribution. A 95% volume contour, for example, typically contains
  95% of the points used to generate the kernel density estimate. Kernel
  parameters were set as follows: scaling factor- 1000000, kernel-
  bivariate normal, single parameter smoothing factor (h)-1000, raster
  cell size-100, and percent volume contours- 50, 90, and 95.


texts: How to export all the layers from Table of Content to shp file?

texts: I have got lots of different layers in the Table of Content panel in ArcMap. I want to export them all in shp format in one go. How can I do that? So far I have been doing it one by one by right clicking on the layer-->data-->export as shp. and it is taking longer time. I want to do them all together.


texts: When selecting by attribute, unique values display only dots (...)

texts: I'm using ArcMap 10.3. I haven't found any information as to why this may be an issue or the cause. Does anyone know why dots show up in the unique value box when I have values in the field in the attribute table?


texts:
I have a marker called chicagoMarker and each time the user zooms in closer I would like to alert them if they are getting closer to chicago. I am not sure if I have a problem with my event listener in this instance.

if (score==0) {
    chicagoMarker.setMap(map)
    google.maps.event.addListener(chicagoMarker,'bounds_changed',function(){
        zoom = map.getZoom();
        if (bounds.contains(chicagoMarker)){
            window.alert("Good Job!");
            }
        })

}


texts: Google Maps layer copyright popup every time map updated / user input

texts:
I have used a Google Maps hybrid layer in Openlayers in my HTML5 iOS app for months now and it worked perfectly. Just recently they changed something because the 'Map data...' line at the bottom looks different (light grey background) and whenever you touch the screen or the map position is updated via navigator.geolocation.getCurrentPosition() the map data copyright box pops up.

Has anyone else experienced this? If so I'd really appreciate any info on any solutions/workarounds you've found.


texts: Where does ArcGIS persist shared AddIn folder configuration

texts:
We want to rollout our customized ArcGIS AddIns enterprise wide. Therefore the esriAddinfiles are stored on a shared folder on the network. ArcGIS does automatically install those (and any updates) if you configure the shared folder in the AddIn Manager. This is described in http://resources.arcgis.com/en/help/arcobjects-net/conceptualhelp/#/Sharing_add_ins/0001000009tm000000/ under "share AddIns on a network"

I wonder if it is possible to configure this folder centralized, like a registry key where it is to set. But I already put a registry tracker on ArcMap to find out where this could be done and found .... nothing. It seems that this is persisted elsewhere.

Anybody has an idea, where I can configure the ArcGIS AddIn shared folder without letting my enterprise end users do it for me (hey, at least I'm sysadmin).


texts:
I'm trying to get python script to open QGIS and display a project file from a windows batch file, but it QGIS never opens. What am I doing wrong?

batch file contains (I've used this with another script which works OK):

SET OSGEO4W_ROOT=C:\filemanagement\Python_QGIS
SET QGISNAME=qgis-ltr
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
set QGIS_PREFIX_PATH=%QGIS%
REM Python Setup
set PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin;%PATH%
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
set PYTHONPATH=%QGIS%\python;%PYTHONPATH%
set GDAL_DATA=%OSGEO4W_ROOT%\share\gdal
"%OSGEO4W_ROOT%\bin\python" "C:\filemanagement\Python_Scripts\QGIS_open_project\open_project.py"

python script below:

from PyQt4 import *
#from PyQt4.QtGui import *
from qgis.utils import *
from qgis.core import *
#from PyQt4.QtXml import *
from qgis.gui import *
import math
import os

from qgis.core import QgsProject
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsApplication

inqgis = iface is not None
#if inqgis == 0:
path="C:/filemanagement/Python_QGIS/apps/qgis-ltr"
qgs = QgsApplication(sys.argv, True)
qgs.setPrefixPath(path,True)
qgs.initQgis()

# Get the project instance
project = QgsProject.instance()
canvas=QgsMapCanvas()

bridge = QgsLayerTreeMapCanvasBridge(QgsProject.instance().layerTreeRoot(), canvas)
project.read(QFileInfo('c:/data/qgis_tests/test.qgs'))
#print project.fileName()


texts:
I have 7 shapefiles, each containing a set of points.

I want to run a 'cumulative' viewshed calculation, which shows not how many points I can see in total, but how many of the 7 shapefiles can be seen simultaneously across the span of my DEM.

Essentially I want to be able to highlight which areas all 7 shapefiles are visible, and where you might only be able to see 1 or 2, regardless of the total number of points contained in the shapefiles.

I looked at this related question, but as far as I can understand it, the cumulative option in the Viewshed tool relates to the cumulative viewshed of all the points within the one shapefile - there is no option to input more than one file. However, I wondered if the answer may lie in running individual ZTVs for each shapefile, and performing some sort of raster calculation with the results?

Can anyone help me out with some advice?

I'm running Las Palmas 2.18.15 on Mac Sierra. I'm a QGIS beginner (no Python experience.)

EDIT: I could merge them, but if I understand correctly, this would result in a viewshed which highlights visibility of ANY of the points, regardless of which shapefile the points originally came from? I want to keep the points grouped, and calculate how many groups of points I can see from any single position - it doesn't matter which points or groups are visible, just if one group, two groups, three groups, four groups, etc. are visible, which is why I have created each group of points in a separate shapefile.


texts: Will hundreds of attributes cause performance problems in an ArcGIS Server layer?

texts:
This relates to ArcGIS Server 10.0 SP1, with data served from a file geodatabase polygon featureclass.

1) For performance reasons, is it a bad idea to have a featureclass with hundreds of attributes?

2) Is it a bad idea to create an ArcGIS Server map service (containing hundreds of attributes) from this featureclass?

3) Is there a performance implication if I build a featureLayer in the JS API, and only specify the attributes I currently need?

eg, the map service might contain 500 attributes, but I define my feature layer using

featureLayer.fields = [x, y, z]

Thanks for any advice, and please let me know if you need any further details.


texts: How to populate one column of a table by extracting the match of a Regular Expression?

texts:
I Have a CSV table where one of the columns has names of the post offices with the type given in brackets like this:

Sri venkateswara puram (Sub Office)

I want to extract the post office type, by using regular expressions, and I can use  (?<=\().+?(?=\)) as the regular expression, but I can't figure out how I can copy the matching string to another column in Qgis's Field Calculator.

How do I extract the match of a Regular Expression in Qgis's Field Calculator, either using the inbuilt functions or custom functions?


texts: Unable to update some of the attributes using PyQGIS


texts: What causes difference between two sources?

texts:
I have downloaded ESA global land cover tiff file and opened it with QGIS 3.6. To measure Africa separetly I used polygon file of this continent as a mask. I got the result that total area of Africa is 29882300 km2 using r.report plugin.
Africa is 30370000 km2 reported in Wikipedia. Why did I get this distinction and how can I reduce it?
Difference is 487700 km2. I think this divergence in values is unacceptable.


texts: JavaScript Api geo processing - internal server error 500

texts:
I am trying to run a geo-processing task using Java Script API for ArcGIS. I have created a simple buffer task, published using ArcGIS server (synchronous) and used the same (almost) code as the view shed example . The view shed runs fine, but it fails with my REST service with the following error.

  Server Error - Object reference not set to an instance of an object.
  Code: 500

.

I have browsed a lot and did some changes with respect but nothing seems to work.
Has anyone experienced this?


texts:
I am using ArcPy to save a TIFF file in a geodatabase file with this code:

arcpy.RasterToGeodatabase_conversion(Input_Rasters=r'D:\Projects\GEBCO_2019_Geotiff.tif', Output_Geodatabase=r'D:\Projects\Example.gdb')

This will create a layer named GEBCO_2019_Geotiff.

What I want is to set by myself the layer name, is there a way?


texts: Undo edits based on sequence rather then layer selected in QGIS

texts: I am used to undo feature of ArcMap that undo everything in sequence irrespective of the layer which is selected. In QGIS it only undo the edits made in the layer selected. Is there a way to change setting so as to undo the edits in sequence irrespective of the layer selected?


texts: Identifying shared polygon boundaries in ArcGIS Desktop

texts:
I am aware you can dissolve boundaries (standard feature analysis tool) but is there a way you can calculate polygons that share a boundary?
I have attached picture with an example of the polygon data I have.


texts: Need advise on map site design (When to wms/kml)

texts:
I'm building a web system from scratch (based on openlayers). All layers in map will be my own (raster and vectors). I want to show polygons/points on map and when user click on them he gets info or list of links of files to download (and need to have the ability using some kind of back office to update new vector layers). I can set up wms server (geoserver or ms4w).
I'm not sure If all layers should be shown as wms or the vectors (which are created in qgis) should be kml files (or something else) and how/where to save the vectors information , in db or on the shp file.

Any recommendation on best practice or guidance for such task will be welcome


texts: Creating input data for LRS plugin from Excel table without coordinates?

texts:
I'm working with QGIS 2.2.0. I want to add a shapefile from an excel file that has informations I want to use in the LRS plugin. The excel file has no geographic information. Every time I try to save it as a shapefile, the program saves in a diferent format that I can not use in the LRS plugin. I import the file as a text layer, with no geographic information, then I save it as shapefile, but that doesn't happen.

To use the LRS plugin, I have defined a road, made by a single line with a length. My Excel file has points with informations about specific kilometric points of that road, but don't have the coordenates, they have the length of the road matching that point.

Do I have to add the respective coordinates?


texts: Fixing OpenLayers Popup position in RTL pages

texts:
I'm using OpenLayers V6 popup but It has wired behavior in Right-To-Left page, When I click on my feature my popup shift map to let and it will be appear in right page! But in LTR page It work fine and when I click on feature It will be appear exactly top of feature, How I can fix this OpenLayers popup in RTL pages, I did change position value to fix in ol-ext(openalyers switcher css) and ol.css but It didn't work.

as demonstrated in image my map is left side but my popup show right side!
This is my code, I inspired from OpenLayers document:

var container = document.getElementById('popup');
var content = document.getElementById('popup-content');
var closer = document.getElementById('popup-closer');

map.on('singleclick', function (evt) {
    //document.getElementById('popup-content').innerHTML ='';
    $('#popup-content').html = '';

var queryLayers =''
         var coordinate = evt.coordinate
  //$('#popup-content').html =''
   var viewResolution = map.getView().getResolution();

      wmsLayers.forEach(function(layer, i){
           if (layer.getVisible()) {
           const value = layer.getSource();
           queryLayers = `${value["params_"].LAYERS}`;

const url = value["getFeatureInfoUrl"](
    evt.coordinate,
    viewResolution,
     'EPSG:3857',
   {'INFO_FORMAT': 'text/html','QUERY_LAYERS':`${queryLayers}`},
);

if (url) {
    //headers.set('Authorization', 'Basic ' + window.btoa('username' + ":" + 'password'));
    fetch(url, {
        method: 'GET', // *GET, POST, PUT, DELETE, etc.
        mode: 'cors', // no-cors, *cors, same-origin
        cache: 'no-cache', // *default, no-cache, reload, force-cache, only-if-cached
        credentials: 'same-origin', // include, *same-origin, omit
        headers: new Headers({
        'Content-Type': 'text/html',
        'Authorization': 'Basic ' +btoa(`${geoUserName}:${geoPass}`),

    'Content-Type': 'application/x-www-form-urlencoded',
    }),
    redirect: 'follow', // manual, *follow, error
    referrerPolicy: 'no-referrer', // no-referrer, *no-referrer-when-downgrade, origin, origin-when-cross-origin, same-origin, strict-origin, strict-origin-when-cross-origin, unsafe-url
    }).then(function (response) {
         //console.log(response)

         return response.text()

      }).then(function(html) {

   $('#popup-content').html(html)

          });

When I change position in DevTools to static popup will show correctly above selected feature.


texts: Split a shapefile into multiple shapefiles based on a particular Field

texts:
I have a shapefile where there are 3 fields OBJECTID, NAME , LOCALITY. I want to create new shapefile that contains the geometry, projection etc. of the parent shapefile and according to the locality;

OBJECTID    NAME    LOCALITY
1    JOHN    UK
2    ANITA   USA
3    CAL     USA
4    DAVE    UK

This shapefile is read and two seperate shape files are created with the properties
of the parent shape files

OBJECTID    NAME    LOCALITY
1    JOHN    UK
2    DAVE    UK
OBJECTID    NAME    LOCALITY
1    ANITA   USA
2    CAL     USA

Please guide me how to do it. I am new to Python


texts: WMS GetCapabilities and variable substitution with MapServer

texts:
I have created a layer with a variable substitution to request wether or not a bird species has been observed in a town. It works ok.
But if I make a WMS GetCapabilities I only see one layer as expected.
I have more than hundred species so if I want to have the GetCapabilities request to show one layer for each species do I need to generate a huge mapfile?
Does anyone can suggest some tips/trick to make this work easier?


texts: Error after executing hex grid script in QGIS 2.6 for Mac

texts:
I'm trying to load a script with QGIS, hex grid from layers bound, with a population shapefile (popp.shp in egis sample data folder). After selecting the layer and executing the script it gives me an error.

Where is the problem?

I attach a couple of screenshots.


texts: Dynamically updating path weights with route in ModelBuilder?

texts:
I am trying to force the routing to reuse as many path parts as possible as described in initial question. The submodel for it is as follows:

What it supposed to do is as follows:

The path from the first iteration is precalculated with normal "Route" in network analyst and given to this model as initial value for "previously chosen route"
Then the iterator goes through the rest of the points and calculates the shortest path under scaled cost (0.1) for the already chosen route segments
The feedback is supposed to dynamically update the chosen route list

The add scaled cost is modeled as follows:

There I have also chosen the option to Append to existing to make the list of previously chosen routes.

The problem is that this dynamic update does not work and I do not know, why.
The model starts showing that there is not enough input parameters and at the same time Copy features does not create any new layers, so practically only the initial input route stays and the rest of output is not created. The iteration and merging is working fine.

The model works, when I remove the feedback link, but this is not what I want.


texts: Placing layer name beside their symbols in legend using representations of ArcGIS for Desktop?

texts:
I want to place the layer name of my feature class (styled with a representation) beside (right hand side) its symbol in the legend of my map in ArcMap 10.1.

I know that the format can be changed via: legend > properties > theme > choose my layer and click "style" and then changing the legend elements' style from the list. However, if I use feature classes styled with representations I do not know how to do this. Using shapes it's no problem because I can adjust the description of the layer in the table of contents. However, this is not possible (or I don't know how to) when using representations.

These pictures visualize my question:

Standard legend style for Representations:

Possible by changing the legend's style:

What I want to achieve:

Any hints or suggestions on how to solve my problem? Maybe there is a way of creating customized styles?


texts: Create a point from a lat/long list using GDAL

texts:
I have a shapefile with a bunch of wildfire perimeter data, and a dataframe with lat and long columns about the lat and long of fire origin points. The information is not connected except by date, so for each fire perimeter, I am trying to check which origin points occurred on that day and fall inside the shape. To do so I open the csv files as dataframes o and p, and then I open the shapefile with

driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open(r'C:\Users\Main\Desktop\locs\here\firep18_1.shp')

# Get Projection from layer
layer = shp.GetLayer()
spatialRef = layer.GetSpatialRef()
print (spatialRef)

# Get Shapefile Fields and Types
layerDefinition = layer.GetLayerDefn()
print(layerDefinition)
inFeature = layer.GetNextFeature()

Then I run the following code.

while inFeature:
    shape = inFeature.GetGeometryRef()
    print(shape)
    break
    """ get the date attribute for the input feature """
    date = inFeature.GetField('ALARM_DATE')
    date = pd.to_datetime(date)

    if date:
        if int(inFeature.GetField('CAUSE')) == 14:
            incidents.append(inFeature.GetField('INC_NUM'))
            temp = o[o['ALARM_DATE'] == date]

            lat, long = list(temp.LAT83), list(temp.LON83)
            for i in range(len(lat)):
                x = lat[i]
                y = long[i]
                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(x,y)
                print(x, y, point, sep='\t')
                if shape.Contains(point):
                    print("something")
                    break

    # destroy the input feature and get a new one
    inFeature = None
    inFeature = layer.GetNextFeature()

However, the point is not being created correctly on lines

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(x,y)

I have tried creating a small contrived example here

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(-121.4347, 39.8134)

poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddPoint(-122, 40)
poly.AddPoint(-120,40)
poly.AddPoint(-120,39)
poly.AddPoint(-122,39)

print(type(poly))
print(poly.IsEmpty())
print(poly.Contains(point))

which also does not work, even though I am sure the point should be within the polygon. Could anyone help me understand what AddPoint takes as arguments, because I think giving it lat long in the form (-121.4347, 39.8134) is incorrect, but I could not find useful examples or documentation.

Edit: When I run shape = inFeature.GetGeometryRef(); print(shape) I get
POLYGON ((138036.905699998 -402646.3633,138086.986 -402651.0541...138036.905699998 -402646.3633)) But I am not sure how they get these values. At first I thought it was just lat and long * 1000, but the multiplied value is different for long and lat.


texts: Is there a 9 digit Zip Code (Example- 75001-6664, Zip + 4) data set that I can download? Please advise.


texts:
I wanted to test the r.viewshed module from grass via the Toolbox in QGis. But no File is generated. This is the log output:

Algorithm r.viewshed - Computes the viewshed of a point on an elevation raster map. starting...
g.proj -c proj4="+proj=longlat +datum=WGS84 +no_defs"
r.external input="/home/matt/Desktop/slopes/N34E069.hgt" band=1 output=tmp14931314968715 --overwrite -o
g.region n=35.0004166667 s=33.9995833333 e=70.0004166667 w=68.9995833333 res=0.000833333333333
r.viewshed input="tmp14931314968715" coordinates=69.6949645541,34.744452777 observer_elevation="2.00" target_elevation="0.0" max_distance="5000" refraction_coeff="0.14286" memory="500" -b output=output1415b59dffed4d07a82eb7ffa7b96a65 --overwrite
g.region raster=output1415b59dffed4d07a82eb7ffa7b96a65
r.out.gdal --overwrite -c createopt="TFW=YES,COMPRESS=LZW" input=output1415b59dffed4d07a82eb7ffa7b96a65 output="/home/matt/.qgis2/processing/outputs/viewshed_test.tif"
Starting GRASS GIS...
Executing ...
Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default
Projection information updated
WARNING: Over-riding projection check
Reading band 1 of 1...
r.external complete. Link to raster map created.
Computing events...
Aborted
ERROR: Raster map not found
ERROR: Raster map or group not found
Execution of finished.
Cleaning up temporary files...
Starting GRASS GIS...
Executing ...
ERROR: Raster map not found
ERROR: Raster map or group not found
Execution of finished.
Cleaning up temporary files...
Converting outputs
Loading resulting layers

The following layers were not correctly generated.
Intervisibility
You can check the log messages to find more information about the execution of the algorithm

It says that the raster map was not found, however, the file exists at the specified location. I am using QGis 2.18.6

What am I doing wrong?

EDIT:
This is the Dialogbox for the r.viewshed:

I've also tried with the python console, but there I do not get any log messages (no output file either):

import processing
from PyQt4.QtCore import QFileInfo

rasterPath="/home/matt/Desktop/slopes/N34E069.hgt"
fileInfo=QFileInfo(rasterPath)
baseName=fileInfo.baseName()
rLayer = QgsRasterLayer(fileInfo.filePath(), baseName)

# Define extent
extent=rLayer.extent()
xmin=extent.xMinimum()
xmax=extent.xMaximum()
ymin=extent.yMinimum()
ymax=extent.yMaximum()

extStr="%f,%f,%f,%f" % (xmin, xmax, ymin, ymax)

# define viewpoint and outputDir
coordStr = '%d,%d' % (34.841101,69.679670)
outputViewshed = "/home/matt/Desktop/slopes/viewshed_test.tif"
processing.runalg("grass7:r.viewshed", rLayer, coordStr, '0', '0', '800', 0.14286, 500, False, False, False, False, extStr, 0, outputViewshed)


texts: I have a raster of forest and non-forest data. I am trying to analyse it using LeCos on QGIS. The pixel size is equivalent to 1sq.m and there class information ranges from 1-503. When I load it into LeCos it doesnt generate any information and QGIS stops functioning. Could someone guide me as to how I can work with LeCos and prevent the software from crashing?


texts: Should ArcGIS change the path of a gdb to the server data store when publishing a GP service?

texts:
I've published a GP service to ArcGIS server and it fails when I try to run it. The workspace is a file GDB on the local network, and the script contains a UNC path to it. There is a copy of the data on the server, and this GDB is referenced in the data store. See screenshot below.

Now I'm trying to call this service using the Javascript API and the service fails with "ERROR 000732: Input Raster: Dataset AE010S010_020M6004P01M010 does not exist or is not supported."

I've found the python script on the server and the workspace path is still set to the UNC path in the original script. The problem is that when running the GP service, the user is local to the server and does not have access to domain resources, so cannot see the network GDB.

The question is, should ArcGIS change the workspace path in the script to D:\data? I thought that was the whole point of adding a data store entry to the server.


texts: Polyline will not turn into a feature layer via MakeFeatureLayer_management (arcpy)

texts:
I've got a python script that is connecting a sequence of points to make lines (script is based on answers to my previous question). To attempt to make subsequent processing more efficient, I would like to do analysis/processing to each polyline as the points are connected, rather than storing them cumulatively in a shapefile (or feature class) and then reading them back into the script later. So the revised code (just the functional looping part, for the sake of brevity) is basically this:

for start in start_cursor:
    startFeature = start.getValue(shapefieldname)
    pnt1 = startFeature.getPart()
    for end in end_cursor:
        endFeature = end.getValue(shapefieldname)
        pnt2 = endFeature.getPart()
        point.X = pnt1.X
        point.Y = pnt1.Y
        array.add(point)
        point.X = pnt2.X
        point.Y = pnt2.Y
        array.add(point)
        polyline = arcpy.Polyline(array)
        arcpy.MakeFeatureLayer_management(polyline, "sightline_lyr")

However, MakeFeatureLayer_management doesn't seem to exactly be working. Any subsequent functions I want to apply (e.g., FeatureToRaster_conversion) to that new feature layer that was supposed to be created from that polyline has no data. If I add the new feature layer to the data frame while the script is running, it shows up in the table of contents in the appropriate data frame, but with the little red exclamation mark that indicates ArcMap can't find the data for that feature layer.

Is an arcpy.Polyline able to be made into a (temporary) feature layer? If not, how can I make a feature layer such that I can do more spatial analysis in further steps? (Short of "make a shapefile or feature class, then read it back in," preferably...)


texts: QGIS raster from interpolation not scaled correctly

texts:
I am using the interpolation plugin in QGIS to interpolate values within a boundary.

I start with two layers; 1) my points and 2) a polygon for the area of interest. Then I go to Raster>Interpolation>Interpolation and add my vector points with the corresponding attribute as an input. I use IDW (TIN does not seem to work for me), set the output file, add result to project, and hit OK.

This creates a raster but it is not the same scale as the point data. Specifically the raster is shorter and narrow (proportionally) than the point data, so the points extend beyond the edges of the raster but the values of those points are represented in the raster. It really just looks like the raster is a little smaller than it should be.

Has anyone experienced this problem? I am not sure what causes it or how to fix it.

Edit:


texts: QGIS Crashing in Ubuntu / 'Save As' browsing

texts:
When trying to create a new shapefile in QGIS on Ubuntu, I get a crash whenever I browse to a folder to save the file.

Has anyone else experienced this?


texts: How to make a virtual field function use the same Ellipsoid settings as the normal field calculator?

texts:
I have a new shapefile layer in CRS 3857, in a project that has the same CRS.
Under project properties, the measure tool is on ellipsoid WGS 84.

When i use the field calculator to create a new field for the layer and use the function "$area", the fields get filled with the correct values.

Now when i use the field calculator to add a virtual field with the function "$area", the values in the field represent the area value of the features calculated without ellipsoid (None/Planimetric).

Is there a way to make the virtual field function use the same Ellipsoid setting as the normal field calculator?


texts:
I would like to know what kind of computation I could perform on a terrain database.
I'm making terrain for VBS2 from raw data (shp files and DTM Elevation files) and then I'd like to be able to check correlation between these two terrain. I want to check that VBS2 created terrain is ''enough'' correlated with the previous data.
I planned to check the shp files and DTM files of course.
But I was wondering if there is mean, median or coefficient that I could calculate on the terrain. Moreover, should I grid the terrain and analyze each square?
Thanks


texts: GeoPandas: Find nearest point in other dataframe

texts:
I've got 2 geodataframes:

import geopandas as gpd
from shapely.geometry import Point
gpd1 = gpd.GeoDataFrame([['John',1,Point(1,1)],['Smith',1,Point(2,2)],['Soap',1,Point(0,2)]],columns=['Name','ID','geometry'])
gpd2 = gpd.GeoDataFrame([['Work',Point(0,1.1)],['Shops',Point(2.5,2)],['Home',Point(1,1.1)]],columns=['Place','geometry'])

and I want to find the name of the nearest point in gpd2 for each row in gpd1:

desired_output =

    Name  ID     geometry  Nearest
0   John   1  POINT (1 1)     Home
1  Smith   1  POINT (2 2)    Shops
2   Soap   1  POINT (0 2)     Work

I've been trying to get this working using a lambda function:

gpd1['Nearest'] = gpd1.apply(lambda row: min_dist(row.geometry,gpd2)['Place'] , axis=1)

with

def min_dist(point, gpd2):

    geoseries = some_function()
    return geoseries


texts: Increase the dimensions of each cell in a raster without any kind of resampling?

texts:
I have a DEM interpolated from a TLS point cloud at 1 mm cell resolution.  This resolution is too fine for use in a 2D surface inundation model with the required timestep, without the model taking far too long to run.  I would like to know if there is any way to proportionally increase the size of each cell and its z value in the raster?

For example, if a cell is currently 0.001 x 0.001 m, is it possible to multiply these dimensions and the z value for each cell by 1000, resulting in a raster that preserves the morphology of the DEM but is at 1 m resolution?

Note also that each DEM is referenced to a local coordinate system set up in a lab.

I'm working in Arc 10.4.1 desktop (advanced licence), but happy to try other open source GIS options.


texts:
Question from my father who is looking to explore alternatives to Manifold:

"I am an independent environment consultant who is a heavy user of GIS for spatial data crunching, geo-referencing self-digitised maps, terrain modelling using LIDAR and analysing GPS field data. I have a workstation 8 core Intel i7 3.4 GHz, 16 Gb of RAM, NVIDIA GeForce GTX 570 running 64 bit Windows 7 ultimate and dedicated HDs for system and cache files. For many years Manifold GIS has been my software of choice due to intuitive GUI, wide range of functions, programmability, natively 64 bit and low price. However, the current version (Manifold 8) is now 4 years old and is now very slow for many tasks. See here for a recent vector processing comparison where Manifold comes out 40x slower than ArcGIS. While Manifold has an impressive suite of surface rendering and analysis functions it is so slow in processing LIDAR data that is almost unusable for anything other than very small data sets. Manifold is CUDA enabled for a limited set of terrain processing functions (see the manual > “index” > “N” > “NVIDIA”) but not the actual rendering. One reason why the next version Manifold 9 is taking so much time may be because the big chunks of software are being coded for CUDA which everyone acknowledges is difficult. Manifold’s loyal user base has been promised the massively faster and improved Manifold 9, but there has been no sign of this new release in 2 years. I, like many Manifold users, are beginning to despair of any development of the software. To protect my business and keep ahead of the game I am being forced to look for other GIS solutions.

An alternative to Manifold would have to have the following:

Run natively in 64 bit and so address large amounts of RAM (ie <
4Gb and typically 16Gb and above).
Able to use the GPU for image/surface intensive processing eg CUDA on
Nvidia graphics cards.
Terrain modelling and LIDAR processing capabilities, be able to
‘drape’ vector and raster map components and do 'fly throughs’ for
video presentations.
Be a stable release.

Manifold 8 professional with ‘surface extensions’ costs $590 (£370). None of the options that I have seen so far seem up to the job. These include:

ArcGIS – I have just been quoted  £1495 + £425 per annum to access
upgrades and    another £2,495 for the 3D Analyst extension.
QGIS – Seems only to support 32 bit in the stable 1.8.0 Windows
version, buggy, and lacking in-house digital terrain mapping
features.
R's spatial packages – this option seems to be too fiddly and static for what I want to do.

Looking forward to responses from people who have faced similar dilemmas – what's the best solution out of these options or should I go for something else completely?"


texts: How to add sfcgal to an already postgis enabled database

texts:
I have a postgres 9.1 database with postgis 2.1, and I recently needed to add sfcgal to it, so i rebuilt postgis 2.1, but I can't find out how to update the postgis extension in postgres without dropping all data.
Is there a way to do that ?

ALTER EXTENSION postgis UPDATE;

returns a notice saying i'm already at version 2.1.0

Thanks for any tips


texts: In FME I want to copy all the styles from a dataset but not the geometry or attributes

texts:
I tried reading it in and exposing all the format attributes and then I connected it to the writer that does the last step of a larger process. It copies the formatting across but also the features. I also tried to put the styler in between.

How can I just get all the styles to use without having to use Mapinfo styler and defining all the styles manually?

Does FME have something like https://github.com/NathanW2/MapInfo-to-QGIS-style-generator but to go from one mapinfo tab file to another?

Can it also copy hotlink settings?


texts: How can I perform linear regression on a file geodatabase table in ArcGIS?

texts: I'm trying to use the Ordinary Least Squares tool in ArcGIS 10.2 but it only accepts feature classes as input and I have a file geodatabase table. I realize that I can export the table to Excel or MATLAB and perform the linear regression there, but I am working in ModelBuilder and need to keep everything within the model. Are there any workarounds or conversions that anyone knows about?


texts: Improving quality of tiles generated from Ordnance Survey Raster Data

texts:
I'm trying to generate some map tiles based on the Ordnance Survey Open Data (available from here).  I don't have much experience with this, but I've read a few tutorials online and have managed to produce some tiles.

The process I'm using is:

# Convert to expanded RGBA
for i in *.tif
do
        gdal_translate -of vrt -expand rgba $i $(basename $i .tif).vrt
done
gdal_merge.py -o merged.tif *.vrt
gdal2tiles.py -r METHOD -z 12 -e -w leaflet -t "Ordnance Survey Raster 250k" merged.tif output_folder

I've tried with METHOD as antialias and as lanczos, although I've done smaller experiments with all the other options.

It works and produces tiles, but they look terrible, so I assume I must be doing something wrong.

The first image below is the result of the lanczos generation, with the tiles combined for a small area around Bath (Zoom: 12, X: 2020-2022, Y: 1363-1365 in google format or 2730-2732 in TMS format). With antialias as the method, the image looks almost identical.

By contrast, below is the same area that I downloaded a few years ago from (I think) cz.tileserver.com/osnew, although that doesn't seem to exist anymore, so I can't check whether that's right.  As you can see, the image quality is much, much better than what I'm able to produce.

Can anyone point me in the direction of what I need to be doing differently?  Or is it just a case of the warping being a fundamental problem and maybe the tiles I've seen before were un-warped and hence not quite in the right place geographically?

Edit: as noted in a comment below, it's more important to me that the map looks good than that it's exactly right in terms of projection, so I'd be happy to find a method that generates tiles and puts them in as close a place as possible if that's always going to look better than the warped tiles.  I know what I'm trying to achieve is possible as the lower of the two images above was created from tiles downloaded from a map server.  That map server had achieved what I want to do, I just want to be able to do the same myself from an OS Open Data download.


texts: Exporting identified intersections from mapshaper

texts:
I have a 3000 mile long line that has quite a few gps errors where the line self intersects itself. Mapshaper does the best job of all the tools identifying these intersections on import, QGIS's intersections tool does not identify them.

I can't find anything in the commands or documentation of how or where it stores the points it identifies as these intersections. I would like to export the points to be able to fix the errors. Fixing them in mapshaper causes to much simplification to the line itself in other areas.


texts: I'm new to the field of gis but I was wondering if there exists anything in ArcScene for indoor-navigation. I have to model a building and then implement a wayfinding component for it. Does ArcScene provide this feature for me or should I think about switching to another tool?


texts: Determining whether point is on left or right side of road using ArcPy?

texts:
Is there any way to determine the left or right side of road using arcpy? I was able to find if road is going from north to south,south to north,east to west or west to east.link that helped:-
How to determine right side of the street vs left side?.But still this do not identify the sides on which the point lies.

I am using ArcPy from ArcGIS for Desktop with a Basic level license.


texts: Postgis Clipping two postgis vector layer without changing attribute table

texts:
What I want to do is make the purple layer at the back having the same shape as the transparent layer.

However, I don't want any attributes of the transparent layer to be joined to the purple layer, is there any command I can use in Postgis to achieve this?


texts: QGIS splits line feature at all intersections?

texts:
I have a long GPS track that I've converted to a shape file. I need to divide it into segments. The first time I use the "split features" tool, it makes a cut where I draw my line, but also splits the line into separate features at every intersection.

So suddenly I go from having 1 feature, to having 1000... I don't want this, I want to split my one line feature into two at the point I chose. NOT a every intersection.

Is there a setting somewhere I'm missing, or should I file a bug report?


texts: PostGIS: proper way to model object that follows sections of another object's path?

texts:
First, apologies if this has been asked before. I don't really know how/what to search for, and am struggling to just put this question into words.

I'm trying to model aerial utilities in a PostGIS database. I'm getting stuck on how to model what would seem to be a line following sections of multiple other lines and how to represent that.

Most fiber optic cable needs support, and so metal wire rope ("strand") is attached to the poles and then the fiber is attached to that. Other fiber optic cable is designed to support it's own weight and can be directly attached to poles.

Here is a (not to scale!) map that tries to show some of the "test cases" that I need to be able to model:

One of the use cases we want to be able to support is calculating a location on the map when someone wants to find "X meters from the end of fiber cable Y" as some diagnostic/testing equipment locates faults that way. Fiber cables do have footage (or "meterage" - is that a word) printed on them, so that is data that can be collected.

I've tried to figure this out on my own, and here is a screenshot of what I have as a Database schema so far:

The most basic location related piece seems to be a utility pole, with a location and a few other columns (owner, label/tag, etc).

So normally, a length or piece of metal wire rope is attached to multiple poles, at a (different) height on each pole. So I have a table for pieces of metal wire rope (mostly just an ID column). And I have a table for "strand attachments" that should link that strand to a pole (with a column for the attachment height). I suspect that my strand_attachement table also should have a integer or real to sort the strand attachments for a piece of strand into an order that can be shown as a line (using the location points gathered from the poles).

The place I'm running into trouble is that I have a table with fiber cables. Those cables will follow the stand support lines, BUT they don't necessarily stay attached to only one strand line NOR do they necessarily  follow the entire path of a piece of strand line. How should I model & record this in my database? I would assume part of the answer is another many to many table linking strand lines and fiber cables, but how do I represent and store what sections of the strand lines the fiber cable follows? When a cable follows multiple strand lines in its path, how do I record what order it follows the strand line "sections" in (maybe a column representing meters along the fiber cable where it meets the route of the strand section?)

I'm not new to database design, just GIS database design. I'm familiar with foreign keys, many to many vs many to one, etc.


texts: Exporting attribute table that has many rows in QGIS

texts:
I checked previous questions that could be related with this question but none of them gave me an answer.

I need to know how to export attribute table that has many rows inside. I have over 2.5 million rows in the attribute table and I cannot export it to excel or as .txt format.

After exporting that data, I will work on it with Microsoft Access Software.

If there is a way to do it, please someone explain it.


texts: PostGIS make line from shortest distance from points to polygons

texts:
I have two tables. An exploded points parcels layer

create table a as
select id, (ST_Dumppoints(geom)).geom as geom FROM tableblah;

and another parcels table call it b which is a polygon

I wrote a a query that calculates the shortest distance from the exploded parcels/as points to layer b

select id,min(st_distance(a.geom,b.geom)) as min_distance
    from a join b on a.id = b.id group by a.id

Now I want to be able to create that minimum line from the points to the parcel. I am unsure how to format this part. I know I probably need st_makeline.

should I convert the layer b parcels to points as well?

UPDATE

I converted layer b too a points also

create table b_points as
select b.id,(ST_Dumppoints(geom)).geom as geom from b;

then i ran

select distinct on(a.id) a.id, st_makeline(a.geom,b_points.geom) as geom
from a join b_points on a.id = b_points.id order by
a.id,st_distance(a.geom,b.geom)

this gives me this which is somewhat correct, but the parcel I converted to points only produced 4 points...

ANSWERED:

ended up using st_shortestline too replace my use of st_makeline. I converted both my polygon parcels table to lines and ran it.as you can see below the results from using st_shortestline is more accurate than st_makeline in regards to finding the shortest path from the parcels boundaries


texts:
I am trying to do something that is very simple, but am having significant problems accomplishing this. I have two types of datasets. One is a set of shapefiles for each state in the US that will create the geographical boundaries of census designated places. The other is a file that contains the latitude and longitude coordinates of certain store locations. It comes in standard Windows formats (csv, excel etc.) What I would like to do is to calculate the distance from the store to the nearest place border. However, I am having significant trouble with the store locations file:

Problem 1: Basic test of data accuracy, having the stores show up on a map.

The data originally was in a comma delimited file. When I followed the standard instructions (File -> Add Data -> Add XY Data…) to input the X,Y coordinates, they never show up on the regional map that I have selected.
They do show up when I choose a different template ( Conterminous United States as a template). However, they show up in the wrong location and as a single point.
I checked the accuracy of the coordinated by using the “Go TO XY” tool. This time it takes me to the correct place, so all the point coordinates are accurate and the format is understood.
I noticed that the XY tool changes my coordinates to a different format. For example, the degree symbol as well as NWES are included in the tool. Changing my data to this format still does not make it show up.
I played around with changing the X and Y coordinates around to make sure that this is not a problem
I read that excel is the preferred format for table files in ArcGIS, so I decided to change them to excel file. However, ArcGIS would not recognize the file. From reading around on the forums, I found that I need to download Excel 2007 driver. Now I can use Excel data, but I have the same issue.
-There are a lot of different coordinate systems to choose from, I played around with all the US based ones (I have US data) and have the same issue with all of them.
When I wanted to just visually see the data on a map, I used ArcGIS Explorer, everything worked fine and was super easy. Not sure why this is so difficult in ArcMap10.2.2

Problem 2: Getting the store location data into the Near tool.

The Near tool does not recognize the text/excel file in which the coordinates are in.
There are several options for choosing which type of data to convert it too. However, here I am not really sure which one is correct. I tried layer and shape and neither of them worked. Due to problems in Part 1, I can’t be sure if any issues stem from incorrect file type.

This is an example of what my data looks like:

X   Y
-92.26960754    30.34570885
-92.25918579    30.37543678
-92.26924133    30.33198738

These are the census shapefiles I am using (2010 Census):
https://www.census.gov/geo/maps-data/data/cbf/cbf_place.html

Any advice? I am quite new to ArcGIS. The online instructions seem really simple, but I am just not getting the results I want.


texts:
How do you make an unclassed choropleth in QGIS 2.x? I'm using 2.4 and can't find an option in a layer's Style tab like the continuous color option that seems to have existed in QGIS 1.x.

I'm curious to see ways to work around this if it is indeed no longer an official option.


texts: ArcGIS Runtime SDK for .NET API method to convert shapefile to ArcGIS feature set?

texts: Could anyone please let me know is there any ArcGIS Runtime SDK for .NET API method to convert a shapefile to ArcGIS feature set object?


texts: ArcGIS how to capture attribute data from another layer

texts:
I have an arc layer with points and another with lines. The points exist along the lines (as stations). Both layers have an attribute for LineName. How do I pass the specific line name data to the points layer?
In my specific example image, the LineName is displayed but the points table has empty LineName attributes, which are the ones I would like to populate.


texts:
I've got database with json stored as an attribute. Is there any way to parse this information and visualize in QGIS? Here is the basic concept.
Example: I would like to display label from JSON.val2 on polyline.


texts: Error when installing QAD plugin: Plugin 'qad' could not be loaded due to an error calling its method classFactory()

texts:
I am trying to install the QAD plugin and it appears this error:
No se pudo cargar el complemento 'qad' debido a un error al llamar a su método classFactory()

AttributeError: 'QVariant' object has no attribute 'split'
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.10/apps/qgis-ltr/./python\qgis\utils.py", line 334, in _startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/Magdalena/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qad\__init__.py", line 29, in classFactory
    return Qad(iface)
  File "C:/Users/Magdalena/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qad\qad.py", line 409, in __init__
    userLocaleList = QSettings().value("locale/userLocale").split("_")
AttributeError: 'QVariant' object has no attribute 'split'

Versión de Python: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
Versión de QGIS: 3.10.13-A Coruña A Coruña, 311f3bbade

Someone knows how can I fix it? I have Windows 10 Pro.

I do not know whether it is related, but when I reinitiate QGIS the QAD plugin started to work (at least I can see the tools). However, I have another problem now. It is not possible to use any of the plugin's tools, because it does not recognise the layers that I need to edit. In the command, it says that the SRC must be a reproject SRC. I did it, but the problem persists.


texts: Creating features(point) grid on polygon

texts:
I am working on GIS based desktop application using C#. I am using dotspatial library in this project.

Now I need to create a grid of features on polygon. This grid cell (rectangle) should be 20*20 Meter Square.

I have worked on it and able to create grid but facing issue regarding to cell size. Whenever polygon size changed cell size also reduced. My code.

         // Polygon Width = 2335
        // Polygon Height = 2054

        int RowsCount = 111;
        int ColumnsCount = 111;

        var maxPointX = Polygon.Extent.MaxPointX;
        var minPointX = Polygon.Extent.MinPointX;

        var maxPointY = Polygon.Extent.MaxPointY;
        var minPointY = Polygon.Extent.MinPointY;

        var dXStep = (maxPointX - minPointX) / (ColumnsCount - 1);
        var dYStep = (maxPointY - minPointY) / (RowsCount - 1);

        var gridColumnsPoints = new double[1000000];
        var gridRowPoints = new double[1000000];

        //   Calculate the coordinates of grid
        var nextPointX = minPointX;

        for (int i = 1; i <= ColumnsCount; i++)
        {
            gridColumnsPoints[i - 1] = nextPointX;
            nextPointX = nextPointX + dXStep;
        }

        var nextPointY = minPointY;

        for (int i = 1; i <= RowsCount; i++)
        {
            gridRowPoints[i - 1] = nextPointY;
            nextPointY = nextPointY + dYStep;
        }

Output

Now when I tried this code on small size of Polygon then grid cell size also decreased.

I know my approach is not correct, So I have searched on it and got some tools. like

Creating spatially projected polygon grid with ArcMap?

But I want to create it in C# and unable to found any algorithm or any other helping material.


texts: How to extract an adjacency matrix from OSM?

texts: I would like to use graph algorithms from the Boost graph library on OSM data. What are my choices to use Openstreetmaps data from within a C++ program ? I need at least a flat file with a list of nodes and edges to build the graph using Boost graph. I am not looking for code that computes anything from the data. I'm thinking of something like an XML parser in C++ that I could integrate to my code.


texts: How to export "solid" 3D objects from QGIS to DXF?

texts:
I work in architecture and we often have to create urban contexts for our projects which involves obtaining a building footprint file and then extruding them by hand.

I have attempted to use QGIS (v3.2.2, macOS) to achieve the goal of generating a DXF file with "solid" extrusions of buildings. My process so far:

Import footprint data
Add Z values using the elevation attribute associated with features
Then process resulting layer with GRASS v.extrude using the approximate height attribute associated with features
Use Export Project to DXF to generate the DXF file

When I loaded the DXF into Revit I expected to see "solid" extrusions. By solid I mean, volumes with opaque faces. Instead, I got 3D wireframe representations of the extruded footprint, as illustrated in the following image:

Is there a way to generate "solid" 3D geometry in QGIS and then export it as DXF preserving the "solidity" of the geometry?


texts: Tools for building and exporting high resolution geotiff

texts:
I'm having some really nice raster maps in geotiff format. I would like to add some features to these maps as well as then exporting them as tiles for serving the map from a web server.

Preferably I would want to be able to work with the different geotiff files in a software like QGis or similar. As I want to add some vector components and combine the different TIFF files before I export them to one big geotiff with all features needed.

Though my biggest problem is that in QGis I can't export geotiffs with higher resolutions than the current screen resolution, which makes it rather pointless as I'm working on a quite big map.

So I'm wondering what possible ways there are to do what I have just described. Can I combine QGis with any other tools or how would you do it?


texts: Earth Engine: Splitting a feature into two using another geometry

texts:
I am trying to split a feature into two based on a line geometry, but couldn't find any function to perform that operation. I have tried cutLines, intersection etc. but to no avail.
Here is my simple code.


texts: Leaflet - click on map object will trigger click on another map object

texts:
i have this situation:

I want that clicking on the circle (buffer) around the marker will trigger a click on the marker.

This is where im creating the circle around the marker:

L.geoJSON(buffer).addTo(this.map).on("click", (e)=> marker.fire("click")

But unfortunately it doesnt work. The click on the circle is triggered but not triggering the click on the marker.

Any suggestions?

Edit: this is where im trying to attach the event:

this.drawnItems = L.featureGroup().addTo(this.map)
let currentLayer = this.drawnItems[index];

let bufMapObject = L.geoJSON(buffer).addTo(this.map).on("click", (e) => {
                console.log(currentLayer)
                currentLayer.fire("click")
            });

The output of console.log(currentLayer) is:


texts: Dynamically selecting t layer features implementing canvasMoveEvent of QgsMapToolIdentify class?

texts:
In my QGIS plugin I need to dynamically select features of a layer with the mouse move event.

In Programatically check for mouse click in PyQGIS? the canvasMoveEvent of the QgsMapTool class is implemented to get x,y coordinates, in Getting attribute value from vector layer after mouse click in PyQGIS? no "move event" is explored (but I think it's not important thing) but they recommended QgsMapToolIdentify to retrieve information about features/layer.

I tried both of these 2 strategies but neither works fine.

Here is my code:

In the dialog implemetation (in the method connected to a button clicked):

    selectTool = CustomSelectTool(self.iface.mapCanvas(), self)
    self.iface.mapCanvas().setMapTool(selectTool)

1) QgsMapTool

class CustomSelectTool(QgsMapTool):

    def __init__(self, canvas):
        QgsMapTool.__init__(self, canvas)
        self.canvas = canvas
        self.identifier = QgsMapToolIdentify(canvas)
        self.featureId = None

    def canvasMoveEvent(self, ev):
        ### IT WORKS BUT MOUSE TRACKING IS NOT VERY ACCURATE
        indentifiedFeatures = self.identifier.identify(ev.x(), ev.y(), self.identifier.TopDownAll)
        if not indentifiedFeatures[0].mFeature.id() == self.featureId:
            self.featureId = indentifiedFeatures[0].mFeature.id()
            indentifiedFeatures[0].mLayer.select(self.featureId)

the cursor doesn't seem to be very reactive, sometimes it misses some feature.

2) QgsMapToolIdentify

    class CustomSelectTool(QgsMapToolIdentify):

        def __init__(self, canvas):
            QgsMapToolIdentify.__init__(self, canvas)
            self.canvas = canvas
            self.featureId = None

        def canvasMoveEvent(self, ev):
            ### NOTHING HAPPEN (IT SEEMS TO BE NEVER CALLED)
            indentifiedFeatures = self.identify(ev.x(), ev.y(), self.identifier.TopDownAll)
            if not indentifiedFeatures[0].mFeature.id() == self.featureId:
                self.featureId = indentifiedFeatures[0].mFeature.id()
                indentifiedFeatures[0].mLayer.select(self.featureId)

in this case nothing happen when I move my mouse on layers in the map. If I put a "print" in the canvasMoveEvent method it doesn't print.

In Capturing mouse movement on QgsMapCanvas? they said there are another two ways to follow:

3) xyCoordinates event

4) eventFilter

but I'm not so expert in PyQGIS developing and no simple examples are shown.

How do I solve this problem?


texts: Find all polygons upstream of a set of points?

texts:
Given two datasets 1) polygons for global river basins and 2) points for global hydro power dams; how can I determine the upstream basins for each dam?

Not sure where to start with this, looking for suggestions. I think I might need to find a river data set that can be used to find all branches up-stream of a given point then use that to select or id the basins?

Student I'm helping is working in ArcGIS Desktop, 10.5.1.


texts: Heatmap algorithm to visualise point diversity

texts:
Can anyone suggest an algorithm to generate a heatmap for visualising point diversity? An example application would be for mapping areas of high species diversity. For some species, every single plant has been mapped, resulting in a high point count, but with very little meaning in terms of the diversity of the area. Other areas genuinely have high diversity.

Consider the following input data:

x    y      cat
0.8  8.1    B
1.1  8.9    A
1.6  7.7    C
2.2  8.2    D
7.5  0.9    A
7.5  1.2    A
8.1  1.5    A
8.7  0.3    A
1.9  2.1    B
4.5  7.0    C
3.8  4.0    D
6.6  4.8    A
6.2  2.4    B
2.2  9.1    B
1.7  4.7    C
7.5  7.3    D
9.2  1.2    A

and resulting map:

In the upper left quadrant, there is a highly diverse patch, while in the lower right quadrant, there is an area with high point concentration, but low diversity. Two ways to visualise the diversity could be to use a traditional heatmap, or to count the number of categories represented in each polygon. As the following images show, these approaches have limited use, since the heatmap shows the greatest intensity in the lower right, while the binning approach would look exactly the same if there were only one category (this could be addressed by increasing the size of the polygon bins, but then the result becomes unnecessarily granular).

One approach I thought of to do this would be to prime a traditional heatmap algorithm by the number of points of different categories within a defined radius, and then using that count as the weight for the point when generating the heatmap. However, I think this might be prone to unwanted artefacts, like mutual reinforcement leading to very sharp results. Also, closely mapped points of the same type would continue to show up as high concentrations, just not to the same extent.

Another approach (probably better but more computationally expensive) would be:

Calculate the total number of categories in the dataset
For each pixel in the output image:

For each category:

calculate the distance to the nearest representative point (r) [probably limiting by some radius beyond which influence is negligible]
add a weighting proportional to 1/r2

Are there already algorithms that I'm not aware of to do this, or other ways to visualise diversity?

Edit

Following Tomislav Muic's suggestion, I have calculated the heatmaps for each category, and normalised them using the following formula (QGIS raster calculator):

((heatmap_A@1 >= 1) + (heatmap_A@1 < 1) * heatmap_A@1) +
((heatmap_B@1 >= 1) + (heatmap_B@1 < 1) * heatmap_B@1) +
((heatmap_C@1 >= 1) + (heatmap_C@1 < 1) * heatmap_C@1) +
((heatmap_D@1 >= 1) + (heatmap_D@1 < 1) * heatmap_D@1)

with the following result (comments under his answer):


texts: Counting Subdatasets in GRIB using ArcPy?

texts:
I'm trying to utilize the extractSubdataset_management tool on a large set of GRB files. Each GRB file has a different number of subdatasets.

Is there any way to execute this tool on all subdatasets contained in each GRB rather than specifying the number of bands?

code so far as follows:

import os
import arcpy

path = "C:/CHRIS/testfolder/testgrib"
gribs = os.listdir(path)
grib = gribs

os.chdir(path)
print os.getcwd()

grb = arcpy.env.workspace = "C:/CHRIS/testfolder/testgrib"

for item in gribs:
    print item
    os.mkdir(os.path.join(path, item + '.subdata'))
    folder = path + "/" + item + ".subdata"
    for i in range(0,115):
        numb=str(i)
        try:
            arcpy.ExtractSubDataset_management(item, folder + "/test" + numb + ".tif", numb)
        except Exception as e:
            print e
            break

so inside the "testgrib" folder are 2 .grb files. One of the files has 112 subdatasets, while the second file has 114 subdatasets. If I specify my range from 0-112, I get an output of 0-112, but in doing so it will omit the remaining 2 subdatasets from the second .grb file.


texts:
I'm publishing shape files in geoserver through curl,my problem is that geoserver didn't recognizing the projection which is in ESRI WKT format(layer was publishing).I searched for this problem i found that ESRI wkt format don't have the "EPSG:4326" .
What i have to do to recognize the ESRI wkt format projection?


texts: Why can't I access a table in a personal geodatabase (.mdb) with arcpy after installing 64-bit geoprocessing?

texts:
I wrote and ran a simple script last year - it cycled through some line features, buffered them and then intersected with some other lines - it all worked fine.

I'm now trying to tweak the file to use a larger buffer, and it's not finding the database I am using to manage the loop and store the results.

This code:

rows = arcpy.SearchCursor("C:\\CEDAR_GIS\\SPEEDY3_GPS_DATA\\shapefiles\\Analysis\\Intersections.mdb\\final_to_school_routes_LINES")

results in this error:

Traceback (most recent call last):
  File "C:\CEDAR_GIS\SPEEDY3_GPS_DATA\shapefiles\Analysis\GISbuff-GPSline_IntersectionTS.py", line 29, in <module>
    rows = arcpy.SearchCursor("C:\\CEDAR_GIS\\SPEEDY3_GPS_DATA\\shapefiles\\Analysis\\Intersections.mdb\\final_to_school_routes_LINES") #opens search cursor on Overlays_to_run
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\__init__.py", line 1133, in SearchCursor
    return gp.searchCursor(dataset, where_clause, spatial_reference, fields, sort_fields)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 359, in searchCursor
    self._gp.SearchCursor(*gp_fixargs(args, True)))
IOError: "C:\CEDAR_GIS\SPEEDY3_GPS_DATA\shapefiles\Analysis\Intersections.mdb\final_to_school_routes_LINES" does not exist

But that table is in that database, and when I ran the script last year there was no error. I have since upgraded to ArcGIS 10.1 with SP1, which I think meant an upgrade to python. Has there been a change in the handling of .mdb files, or am I doing something stupid?


texts: removeSelection() influences removeMapLayer() - Message: Layer not registered

texts:
I'm currently writing a plugin in QGIS. Through an action, I want all the selected Features and a specific Layer to be be removed.

def deleteLayer(self):

    #remove Selection
    for layer in layerlist:
        if layer.type() == QgsMapLayer.VectorLayer:
            layer.removeSelection()

    #deletes Layer
    self.mapRegistry.removeMapLayer(self.layer)

Both sections work fine, if run while the other is commented out. Yet, the moment the code is like above it runs as it should... But I get this message from QGIS:

In German: Layer20171201140534272 :: Layer nicht registriert.

English_Translation: Layer20171201140534272 :: Layer is not registered.

Does anyone know what the message means and why it appears?

Edit: Maybe this additional code is important - layerlist does not include self.layer, because it was removed.

if self.layer in layerlist:
        layerlist.remove(self.layer)

And the message also appears when removeMapLayer() is given the layer.id().


texts: What is the difference between .filterMetadata(), ee.Filter.lte() and .sort() in google earth engine?

texts:
I am currently trying to obtain images that have less than 70% cloud cover (for example) and I have come across multiple variations of the same. Could you tell me what is the difference between these three approaches and which approach might be most suitable for obtaining a set of least cloudy images?

Approach 1:

var collection = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 70);

Approach 2:

var collection = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .sort('CLOUDY_PIXEL_PERCENTAGE', false);

Approach 3:

var collection = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 70);


texts: Is there a Python equivalent to QGIS Raster Calculator's "file1.asc" AND "file2.asc"

texts:
I would like to script a calculation equivalent to QGIS Raster Calculator's

"file1.asc" AND "file2.asc".

I have 18 data raster files which I would like to do the equivalent of ADD to each of five other rasters.

Is there an easy way to script (Python?) this that will take me less time than ploughing through it in Qgis?


texts: PostGis vs Oracle Spatial 11g performance

texts:
I am having a performance issue updating a table in Oracle Spatial 11G while converting 63,000,000 points to WKT format. I also made a comparison with Postgres. Both Postgres and Oracle run in two virtual machines hosted by the same server.

Postgres: processors 12, ram 64GB

create table alldata(geom geometry(Point, 4326));
alter table alldata add column wkt4326 varchar(100);
update alldata set wkt4326 = St_AsText(geom) -- 25 minutes

Oracle 11G: processors 12, ram 40GB, pga + sga = 8GB maximum

create table alldata(geom sdo_geometry); -- e.g.: SDO_GEOMETRY(2001, 4326, SDO_POINT_TYPE(aLon, aLat, null), null, null)
alter table alldata add (wkt4326 varchar(100));
-- **Did not complete in a hour...**
update alldatafiltered set wkt4326 = SDO_UTIL.TO_WKTGEOMETRY(geom);

Why is Postgres outperforming Oracle while updating tables? Which causes may affect Oracle spatial performances?

At the moment I am not leveraging any spatial index.


texts: How to add ArcPy Feature Layer to ArcMap automatically?

texts: With Python, how do I get either the MakeXYEventLayer_management or the MakeFeatureLayer_management layer that is temporarily created added to ArcMap automatically after the tool is run, without using arcpy.mapping?


texts: Map function over ImageCollection in Python 2.7 in Earth Engine

texts:
I am trying to map a function over an image collection in python.  I am able to run same process in JavaScript, but in Python I am having issues.

Here is the code in JavaScript, which returns a JSON full of values:

var smap = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture').filterDate('2019-04-25', '2019-05-17')

var geom = ee.Geometry.Point([-94.440367, 42.1068856])

var res = smap.map(function(image) {
  var reduced = image.reduceRegion({geometry: geom,
                                    reducer: ee.Reducer.mean(),
                                    crs: 'EPSG:4326',
                                    scale: 30})
  return image.set('value', reduced);
});

var vals = res.aggregate_array('value');
print(vals);

However, when I run this in python, I get an empty list:

smap = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture').filterDate('2019-04-25', '2019-05-17')

geom = ee.Geometry.Point([-94.440367, 42.1068856])

def eeAlgo(image):
    reduced = image.reduceRegion(geometry=geom, reducer=ee.Reducer.mean(), crs='EPSG:4326', scale=30)
    image.set('value', reduced)

    return image

vals = smap.map(algorithm=eeAlgo)

vals.aggregate_array('value').getInfo()

> []

Is there some difference between running this in Python that I am missing? Or is this a bug that I should file at https://github.com/google/earthengine-api? I am using Python 2.7.


texts:
I have a data for a project and the attributes had been written in Arabic and this what I found, how can I repair it and get the Arabic texts?

I am using ArcGIS 10.4


texts: Is there a way to access ps.map in python?

texts:
i am wondering whether its possible using python interface to load layers and edit them in Print/Cartographic Composer without displaying, and finally export to PDF?

something that can be done via ArcPy (ArcGIS python module).

thanks for any inputs!


texts:
I have a gpkg file. I trying to open it with spatial_lite gui and cast geometry to GML.
select AsGML(geom) from MyTable

but i got null. I tried another functions AsText(), AsWKT(), st_asText()... but anyway i got null. BTW function ST_GeometryType(geom) returns LINESTRING, but other functions don't work.
My final goal is to cast geometry from gpkg to GML on c#. But it's looks impossible because I cannot even select geometry in  spatial_lite gui.
Do I use wrong queries? What can be wrong here?


texts: How to split lines that are larger than 10km?

texts:
I have a shapefile of line segments representing global rivers. The segments are all different lengths. If a segment is longer than 10km, I would like to divide it in to multiple segments approx 10km long + a segment with whatever is leftover, because not all segment lengths are perfectly divisible by 10....Alternatively, I could divide the long segments evenly into multiple segments, as long as they are all under or near 10km in length. I'm not sure which is easier.

Since this is a global data set, manually splitting things using 'split' in editor is impossible.

Any ideas to do this for my whole dataset?

I haven't used python with GIS before, but it might be necessary here.


texts: How does MERIS collect vegetation height data?

texts: In GlobCover 2005 and 2009 project, the sensor used for spectral data collection was MEdium Resolution Imaging Spectrometer (MERIS). It observes the Earth in 15 spectral bands in visible and near infra-red region. But I am not getting which of those spectral bands is used for vegetation height detection. Land cover classes in GlobCover consist of vegetation distinguished on the basis of height. But without LIDAR, I do not think that data such as tree or shrub height can be determined. Then how is the vegetation height measured in GlobCover classification?


texts:
This is a follow-up to my previous question:

Where can I find recent rainforest (change) land cover data?

I have succeeded in finding data that more or less meet my needs but I cannot interpret them:

http://www.vgt.vito.be/

Examining the NDVI layer in QGIS, I see that the pixels are integer-valued, but their meaning is not clear. You will notice that on the (older) data on the GLCF site define what these sorts of numbers mean. Unless I'm missing something, I can't find those instructions to go with the SPOT data.


texts: Run spatial clustering by inputting coordinates to cluster algorithm

texts:
I've been reading about algorithms for spatial clustering and it's easy to get lost since there are dozens of them. An idea that came to my mind is to use an a-spatial cluster algorithm on a dataset where the lat and long variables are also included as inputs of the model. This way, the algorithm uses the spatial proximity between observations as one additional co-variate in the clustering process.

Question: This idea sounds too obvious and simple to me and yet I haven't seen any book nor manual mentioning this approach, or at least explaining why it shouldn't be used. So I would like to pick your brains on this, what pros and mainly what cons you see in this approach ?

Clearly, this is an unorthodox approach and it's somewhat limited because it considers only euclidean distances and it would work better for spatial points or regular polygons on a regular lattice. Apart from that, I'd like to hear if there are any strong reasons not to use this approach. Any ideas?

Reproducible example below using mclust in R

library(mclust)
library(dplyr)
library(spdep)
library(sp)

# load housing dataset of Lucas county OH/USA
  data(house)
  house <- spTransform(house, CRS("+proj=longlat +datum=WGS84"))
  plot(house)

# get only a sample of points for the sake of code speed
  df <- data.frame(house)
  df$id <- 1:nrow(df)
  set.seed(1)
  id_samples <- floor(runif(3000, min=1, max=nrow(df)))
  df <- subset(df, id %in% id_samples)

# get variables to input in the model into a matrix
  df <- df[ , c('long', 'lat', 'yrbuilt', 'avalue')] # get variables lat, long, when the house was built and its value
  df <- as.matrix(df)

# run Mclust model
  d_clust <- Mclust(as.matrix(df), G=1:9)

# the optimal number of clusters is 8
  cat("model-based optimal number of clusters:", dim(d_clust$z)[2] , "\n")

# plot. The 1st chart on the 1st column shows how the clusters are spatially distributed
  plot(d_clust, what="classification")


texts: Determining if Point is in Polygon using PostGIS?

texts:
I want to check if a point is inside a fenced polygon.

INSERT INTO polygon (name, geom) VALUES ('test', ST_GeomFromGeoJSON('{"type": "Polygon", "coordinates": [[0,0], [0,5], [5, 5], [5,0]]}'));

SELECT * FROM polygon WHERE name = 'test' AND ST_Intersects(geom, 'POINT(2 2)');
SELECT * FROM polygon WHERE name = 'test' AND ST_Contains(geom, 'POINT(2 2)');

Ideally, I'd expect it to return the row inserted as it should lie within the area. But I get 0 rows for both queries.

I also tried with much bigger, and more realistic data (using lat/long from a map), but I simply cannot get it working and the docs don't seem to help either?


texts:
I have a custom wms server
the transform for this wms is UTM and i cannot set googlemap layr for background
for example my feature in wms is Iran country but googlemap background is Algeria

this is my code :

var options = {

        //maxExtent: new OpenLayers.Bounds(35.7008418592904, 51.3777884858342, 35.7520557036613, 51.4272609453396),
        //maxExtent:new OpenLayers.Bounds(506968.03125, 3936165,557431.4375, 3965082.5),
        units: "meters",
        format: "image/jpg",
        projection: new OpenLayers.Projection("EPSG:2058"),
        displayProjection: new OpenLayers.Projection("EPSG:900913"),
        numZoomLevels: 50
    };
function init() {
    var capaGoogle = new OpenLayers.Layer.Google(
               "Google Satellite",
               { sphericalMercator: true,projection:"EPSG:900913",
               transparent: true
               }
           );

    map = new OpenLayers.Map('MapDiv', options);
    wms = new OpenLayers.Layer.WMS(
                'WMS Layer Title',
                '../wms.ashx?dis=t',
                { layers: '', version: '1.3.0', srs: "EPSG:2058", transparent: true },
                { singleTile: true, ratio: 1, displayOutsideMaxExtent: true }
            );
     map.addLayers([wms, capaGoogle]);
 var lon = 3943679.05500279;
    var lat = 538414.111109505;

    var zoom = 11.3;

    center = new OpenLayers.LonLat(lat, lon);
    alert(center);
    map.setCenter(center, zoom);

}

i used transform like

center = new OpenLayers.LonLat(lat, lon).transform('EPSG:2058', 'EPSG:900913');

but that is not worked :(


texts: View domain code values from ArcMap attribute table?

texts:
In the attribute table of a feature class in ArcGIS Desktop 10.7.1:
Is there a way to view data as the domain codes, rather than the domain descriptions?


texts: Making graph from network vector in QGIS?

texts:
I have a network of roads and I would like to obtain a graph made by edges, each of them delimited by a origin node and a end node, such as the picture shows.

I would like to obtain a table made of rows like this

  "edge_id, origin_node_id, end_node_id, origin_node_X, origin_node_Y, end_node_X, end_node_Y"

for example, for the a12 edge: "a12, 48.1234, 50.0000, 48.1335, 49.1111", where the origin node 1 coordinates are (48.1234,50.0000) and the end node 2 coordinates are (48.1335,49.1111).
I used the field calculator with these expressions:

  x_at(0),y_at(0),x_at(1),y_at(1)

to get "origin_node_X,origin_node_Y,end_node_X,end_node_Y".

Now I am searching for a way to obtain the id of nodes which obviously are in another layer.

So I think that I could need a field calculator expression or python code to relate nodes' layer with edges' layer. I think that "get_feature" could be useful but at the moment I did not find the right way to use it. Another option is to make a text output file with python.


texts: identifying single overlaping features from a GeoServer layer

texts:
I'm trying to use GeoServer to display features from a PostGis table the url will look similar to this.

Example URL

Localhost:8080/geoserver/cite/wms?service=WMS&version=1.1.0&request=GetMap&layers=cite:Layer1, cite:Layer2& styles=&bbox……..width=512,heaight=430&srs=EPSG:4326&application/openlayers

&FILTER= (<Filter xmlns="http://www.opengis.net/ogc">
               <FeatureId fid=”Layer1.1234”></FeaturId>
               <FeatureId fid=”Layer1.5678”></FeaturId></Filter>)
            (<Filter xmlns="http://www.opengis.net/ogc">
        <FeatureId fid=”Layer2.345”></FeaturId></Filter>)

When I click on a polygon I need to highlight a single polygon if two (or more) polygons overlap i need to be able to step through each polygon in turn showing just it's attributes.

At present I click on a feature and the attributes for the entire layer are shown
in my example I use PostGis I've also tried this using a shapefile and the same thing happens.

Any help would gratefully received
ThanX Holly


texts: GeoServer s3-geotiff module with custom S3 provider fails

texts:
I'd like to use the s3-geotiff module within GeoServer with a custom provider (Minio).
The minio server is running properly and I could reach the tiffs in the bucket from a variety of s3 tools.

Following the instructions of the docs I've created a s3.properties file with endpoint, username and password strings.

I've placed the s3.properties file in the root of Geoserver_Data dir and also in the security folder.
Unfortunately, using the s3 buckets tiff as a datastore fails. Geoserver always wants to reach the bucket on an AWS S3 instance. Using AWS, everything works as expected.

Any ideas?


texts: using qgis to split geometry in postgis: how to add unique id and copy original data in linked tables

texts:
We have set up a postgis database to manage all our old shapefiles.

There is one table that contains the geometry and some basic information.  There are several tables linked to that in a 1:m relation.

We want to clean up the geometry (snapping, splitting, moving, etc.). We are using QGIS as a desktop-application to do that. As long as we do not split a geometry, there is no problem in the database.

While splitting a geometry, you create a new feature in the basic table.  This one need to get a new unique id. All original linked features in the other tables have to be duplicated with that new unique id as a reference.

In a manual way I could do

INSERT INTO weg_bron ("weg_id", "bron_id", "bronnummer1", "bronnummer2")
SELECT <id new segment>, "bron_id", "bronnummer1", "bronnummer2" FROM weg_bron WHERE weg_id = <original id>;

What is the workaround to do this (semi)automaticly while working with the splitting tool in QGIS?

--- some more information about the tables ---

The basic table is called 'weg':

 CREATE TABLE weg
(
  id serial NOT NULL,
  naam text,
  geometrie geometry(LineString,31370) NOT NULL,
  lengte integer,
  deelgemeente_id integer,
  verharding_id character varying(2),
  categorie_id character varying(3),
  juridisch_statuut_id smallint,
  breedte_id smallint,
  toegankelijkheid boolean,
  beschrijving_toegankelijkheid text,
  zichtbaar boolean,
  trage_weg boolean,
  opmerkingen text,
  eigenaar_id integer,
  inventariseerder_id integer,
  inventarisatiedatum date,
  aanmaakdatum timestamp with time zone NOT NULL DEFAULT now(),
  aanpasdatum timestamp with time zone NOT NULL DEFAULT now(),
  CONSTRAINT weg_pkey PRIMARY KEY (id)
  )

One of the linked table is called 'weg_bron':

    CREATE TABLE weg_bron
(
  weg_id integer NOT NULL,
  bron_id integer NOT NULL,
  bronnummer1 character varying(25),
  bronnummer2 character varying(25),
  aanmaakdatum timestamp with time zone NOT NULL DEFAULT now(),
  aanpasdatum timestamp with time zone NOT NULL DEFAULT now(),
  opmerking text,
  CONSTRAINT weg_bron_pkey PRIMARY KEY (weg_id, bron_id),
  CONSTRAINT weg_bron_bron_id_fkey FOREIGN KEY (bron_id)
      REFERENCES bron (id) MATCH SIMPLE
      ON UPDATE NO ACTION ON DELETE CASCADE,
  CONSTRAINT weg_bron_weg_id_fkey FOREIGN KEY (weg_id)
      REFERENCES weg (id) MATCH SIMPLE
      ON UPDATE NO ACTION ON DELETE CASCADE


texts:
I'm using QGIS 3.4.1 and TimeManager 3.0.2.

I've tried every combination that I've seen suggested and still cannot get a working animation.  Started w/ a CSV import, saved to an ESRI shapefile.  All the data was brought in as a String.  I've since created new fields of both Date and Integer and copied over the date data into those columns and tried using each as the Start Time option in TM.  Also, I've formatted the date as YYYY-MM-DD.

Are there any issues w/ the version combination I have for QGIS and TM?


texts: How can I edit an existing thematic shapefile composed of contiguous polygons in QGIS

texts:
I cannot edit an existing thematic shapefile layer composed entirely of contiguous polygons. I created the layer many years ago using ArcGIS Desktop; the themes represent the depth of water in a shallow harbour https://www.wexfordharbour.info/iChart. My license for ArcGIS has expired and I cannot afford the €200 for a home license so I decided to try QGis. At first sight it seemed perfect. I opened a project and imported the shapefile. I imported a satellite map, georeferenced it and now I want to modify the polygons according to the underlying sandbanks and shallows showing on the satellite image.
I spent a while removing all geometry errors; they were minor and mostly around the outer box surrounding the shapefile. Now I want to edit polygons to change their shape. QGis Reshape Features does not work for me, nothing happens. Creating a polygon totally within an existing one (Fill Ring) works fine but if and of the surrounding polygon is included, even one vertex, the create does not complete giving a warning "Could not create ring since the inserted Ring crosses existing rings".
If I create polygons using Add Polygon, the result is empty and I cannot select it using Identify Features button. Though the outline is apparent it appears to be a polyline rather than a polygon.
Finally, I cannot clip the underlying layer with a newly created Part of Polygon using a button or even the Plugin Vector/Clipper/Clip. It just gives a host of errors.
If I empty the whole shapefile and digitise Polygons and Rings where I want them all seems fine except that I cannot create a Ring within a Ring if it shares any of the boundary.
I have 30 years of experience with ArcInfo/ArcView/ArcGIS so I do know the concepts of GIS but I am stumped. How can I do those?


texts: I have a map with a parcel (property) layer with property information including the owners name and address.  I would like to label all of the properties with a key number/letter and reference it to a table where the rest of the property information is shown.


texts:
Having a strange issue consuming WMS services hosted on ArcGIS Server 10.3.

When the service is cached in ESRI's compact format, the request fails with HTTP 400. In exploded format, it works perfectly. I've tracked down the issue to the CRS parameter. For example:

`http://gis.example.org/arcgis/services/ServiceFolder/ServiceName/MapServer/WmsServer?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&CRS=CRS:84&BBOX=...'

will fail based on the CRS=CRS:84 parameter. If I change that to CRS=3857, it succeeds. Again, this only happens with compact format. I thought it might have something to do with the BBOX units, but the parameters for exploded are identical and it works.

The service was published as EPSG:3857. For what it's worth, the server was installed as 10.2.2 and upgraded to 10.3. I haven't tested through a fresh install of 10.3 or through an older version, but suspect it wouldn't matter. Have an open ticket with ESRI but have yet to hear back.

I've also tried editing the service's cfg with listCustomCRS and listSupportedCRS to no avail.


texts:
I would like to collect a table of equivalent or close enough gdal based functions to replace ArcGIS functions (esp in spatial analyst). Are there functions available using gdal to perform the following ArcGIS functions:

Lookup
Zonal statistics
Focal statistics
Nibble
Iterate through raster VAT
Combine
Mosaic

The key is to have a ready reckoner of such functions, so ArcGIS based code can be relatively painlessly changed to one based on gdal.