texts: Select Layer by Location Produces Empty Output in Python

texts:
I am new to using Python scripts to execute ArcGIS tools. In my first script I am attempting to loop through many Select Layer by Location operations on overlapping feature classes. I want to create a new feature class that contains all of the "Block" features that have their center in the "Coverage" feature. All selecting polygon feature classes are in a "Coverage" feature dataset and all input feature classes are in a "Blocks" feature dataset. Both feature datasets are in a File Geodatabase. The selection completes with no unexepected behavior in a batch run of a functionally similar ModelBuilder Model, but the Python script below produces empty feature classes for the selected features.

My code is below:

#Import arcpy module
import os
import arcpy
from arcpy import env

#Define workspace
env.workspace = "C:\\ArcPrj\\NFHL_Data.gdb"

#Local variables
#Feature dataset containing blocks
Block_fd = "C:\\ArcPrj\\NFHL_Data.gdb\\CensusBlocks"
#Feature class containing coverage polygons
Coverage_fd = "C:\\ArcPrj\\NFHL_Data.gdb\\NFHLCoveragePolys"
#Output location
Output_fd = "C:\\ArcPrj\\NFHL_Data.gdb\\NFHLBLocks"

#Create lists from feature datasets to select from
env.workspace = Block_fd
blockList = arcpy.ListFeatureClasses()
env.workspace = Coverage_fd
coverageList = arcpy.ListFeatureClasses()

#Iterate through each feature class in the Block dataset and the Coverage
#dataset and perform operations

#Change workspace to main gdb
env.workspace = "C:\\ArcPrj\\NFHL_Data.gdb"

for i in range(len(blockList)):
    #Change workspace to main gdb
    env.workspace = "C:\\ArcPrj\\NFHL_Data.gdb"

    #Select the ith feature class in the feature dataset
    blocks = blockList[i]
    coverage = coverageList[i]

    #Make feature classes into feature layers
    arcpy.MakeFeatureLayer_management(blocks, "blocks_lyr")
    arcpy.MakeFeatureLayer_management(coverage, "coverage_lyr")

    #Select blocks that have centroids in the coverage polygon
    #(select layer by loc)
    arcpy.SelectLayerByLocation_management("blocks_lyr", "HAVE_THEIR_CENTER_IN", "coverage_lyr", selection_type="NEW_SELECTION")

    env.workspace = Output_fd
    #Copy selection to a feature class
    arcpy.CopyFeatures_management("blocks_lyr", blocks + "_NFHL")

    #Clear previous layers
    arcpy.Delete_management("blocks_lyr")
    arcpy.Delete_management("coverage_lyr")

As a Python beginner I appreciate any help you can give me. Thank you.


texts:
I am trying to export a shapefile from the ArcMap (Made in both 10.2 and 10.3; same problem) environment to be used in R.  Unfortunately, it seems that the projection information is somehow not retained at some point during this process.  The closest question I have been able to find related to this is Exporting point shapefile from ArcGIS for Desktop to R?, however the user was having problems on the R end of things. I believe my issue may be on the Arc side:

I have manually created a shapefile using the following steps:

Import DEM of region.
Manually georeference a scanned historical map using notable landforms
Within ArcCatalog, create a new shapefile
Select preferred XY coordinate system (WGS 1984 Complex UTM Zone 22N), and specify that it will be a polygon
Enable the Editor toolbar
Edit the new shapefile, as a polygon.  Select points around the extent of the baselayer region of interest.

However, when bringing the shapefile into R, it seems that the only info that is retained is the extent of the object. I have used readOGR, readShapePoly, and read.shapefile but all of the functions leave me in the same situation.  It seems that using read.shapefile brings the information in as a list:

> summer.s <- readOGR(dsn=".", layer="Trynagain")
OGR data source with driver: ESRI Shapefile
Source: ".", layer: "Trynagain"
with 1 features
It has 1 fields
> proj4string(summer.s)
[1] NA
> summer.s
class       : SpatialPolygonsDataFrame
features    : 1
extent      : 494033.5, 532556.9, 7431381, 7452759  (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables   : 1
names       : Id
min values  :  0
max values  :  0

No matter the import method I use, R just can't seem to figure out what's going on with the projection information.  Yet, when I use the same files in Arc, it knows where and how to project properly.  Does anyone have any suggestions for what I may be doing wrong?


texts: Is it possible to customize openFeatureForm of PyQGIS?

texts:
I am a beginner at PyQGIS.

I am trying to create a plugin, which will be used to Add, Delete, and Edit the Feature.

It is more customized than which is provided by QGIS (2.2).

I am capturing the Points using RubberBand API and then Add the attributes for the feature using openFeatureForm(). The openFeatureForm() takes default values as NULL.

What I want is to customize this form. Like, populate some values with a Combobox (to allow a user to select the required values). Also for fields like the date, I want to put DateTimePicker() widget in the form.

So Is it possible to customize the form provided by openFeatureForm()?


texts: Calculate visibility with polygonal visual barriers

texts:
I have the raster file for my digitial elevation model as well as a polygon shapefile with polygons representing different kinds of environments (forest, field, meadow, etc), and wouldlike to calculate the surface of the visible area from several points while taking into account both elevation and environment type - for example, forest polygons should be visual barriers as it is not possible to see anything through forests. But there isn't such an option in the Visibility tool.

I already tried to calculate visibility within a donut raster which forest polygons are removed from. The good point is that forests themselves are not counted as visible when calculating the visible area, but the problem is that with this method they are considered transparent (as usual with no data cells...), and the area right behind them is always considered visible, which is the opposite of what I wanted...

Does someone know how to do this ? I use both Python and the user interface, but prefer Python.

I was also considering the creation of a local raster with higher raster elevation values where the forest polygons are, or to assign a very high elevation to all cells with no elevation data (which would make them visual barriers) but I have no idea how to do this.


texts:
I'm making a Leaflet map, and using the Marker Cluster plugin to cluster a Few Thousand Markers.

On of the attributes in my Data is the State that this marker lies in.

I would like that at the initial zoom level, the map show the Markers clustered by the states, i.e. Each state have one marker cluster.

How can I do this?


texts: Problem while applying trigonometric functions in QGIS 3.1

texts:
Using the Trigonomteric functions in the QGIS 3.1 field calculator gives me unexpected output. For example sin(180) produces the output -0.8011526357338304 instead 0.

What is the reason for this?

How do I achieve expected result using degree as input?


texts:
I'm trying to find the TauDEM Model Builder Models toolbox.

I'm teaching myself how to use TauDEM to delineate watersheds. After running through some examples I started putting the modules together in Model Builder. Before I finished, I realized that someone may have already done this and I found exactly what I was building on pages 10-11 on David Tarboton's website (developer of TauDEM). However, the model that is referenced as available to download no longer seems available.

This site has downloads for the scripts and individual tools, but not the full model builder geoprocessing tool. Anyone know if this is available anywhere? If not, I will


texts: Problem of linear records with zero shape length and no geometry in attribute table

texts:
I  am working with a bicycle and pedestrian facility file geodatabase in ArcMap.  I noticed a some of the records in the shapefiles/feature classes have a shape length of zero in the shape length field in the attribute tables.  These records have attributes but no geometry that appears when I try to "zoom to selected feature."  This makes since because the shape length is zero.

Does anyone have an idea of how these records might've come about?  How can you have records in an attribute table with no corresponding geometry?  I want to say in the past i've been in an editing session adding new features and seen ArcMap created blank features that have no corresponding geometry when I attempt to "zoom to selected feature."


texts:
I have a series of tabs running down the left hand of my page and on each page I have a leaflet map. Each map is based on the choropleth example but I'm running into problems. Only one of the three maps enables me to do the custom info call to put the attributes into the info panel.

I've quickly worked out that each map has to have its own 'map' ID e.g

...var map = L.map('map1').setView([37.8, -96], 4);...

...var map = L.map('map2').setView([37.8, -96], 4);...

...var map = L.map('map3').setView([37.8, -96], 4);...

so I'm guessing that I would need to make each instance of the info.onAdd function unique? I've played around do this but with no luck.

var info = L.control();
info.onAdd = function (map) {
    this._div = L.DomUtil.create('div', 'info'); // create a div with a class "info"
    this.update();
    return this._div;
};
// method that we will use to update the control based on feature properties passed
info.update = function (props) {
    this._div.innerHTML = '
US Population Density
' +  (props ?
        '' + props.name + '
' + props.density + ' people / mi2'
        : 'Hover over a state');
};
info.addTo(map);


texts: Loading MarkerClusters while the bounding boxs changes position and display

texts:
I'm using Leaflet.markercluster plugin and JavasScript to develop a map where i want to generate clusters according to the screen bounds. I have a service that returns me a list of points depending on screen bounds, with lat and lng properties.

Declaration of variables map and cluster.

/* Single marker cluster layer to hold all clusters */

var markers = L.markerClusterGroup({
  spiderLegPolylineOptions: {weight: 0},
  clockHelpingCircleOptions: {weight: .7, opacity: 1, color: 'black', fillOpacity: 0, dashArray: '10 5'},

  elementsPlacementStrategy: 'spiral',
  helpingCircles: true,

  spiderfyDistanceSurplus: 25,
  spiderfyDistanceMultiplier: 1,

  elementsMultiplier: 1.4,
  firstCircleElements: 8

});

for (var i = 0; i < 10000; i++) {
  var circle = L.circleMarker([Math.random() * 30, Math.random() * 30], {fillOpacity: 0.7, radius: 8, fillColor: 'red',
  color: 'black'});
  markers.addLayer(circle);
}

var sudOest = L.latLng(41.512742, 2.410589);
var nordEst = L.latLng(41.590395, 2.489678);

//la sortida del mapa depen de la capa topònims i concretament d'aixo:      map.fitBounds(toponims);
map = L.map("map", {
    zoom: 14,
    crs: crs25831,
    layers: [baseLayers[Object.keys(baseLayers)[0]]],
    zoomControl: false,
    minZoom: 0,
    maxZoom: 6,
    maxBounds: L.latLngBounds(sudOest, nordEst)
}).setView([40.705008, -73.995581],6);

map.addLayer(markers);

Function called on move where load all coordinates points according to the bounding box

function cargarMegaCapes(zoomActual){
    for (var prop in capesDinamiques) {
      if(map.getZoom()>3){
        if(capesDinamiques[prop][0]){

        removeAreaLayer(prop); //clean the other markercluster

        baseLayers[prop] =afegirCapaServeiWeb(capesDinamiques[prop][1]
        +zoomActual._southWest.lng+'&x2='
        +zoomActual._northEast.lng+'&y1='
        +zoomActual._southWest.lat+'&y2='
        +zoomActual._northEast.lat+''); //service that return a list of points

        addAreaLayer(prop);//crating the new marker clusters according to the bounding box

            }
        }
    }
}

    map.on('moveend', function() {
        var zoomActual = map.getBounds();
        cargarMegaCapes(zoomActual);

    });

Services example (JSON):

e {_layers: {…}, _initHooksCalled: true}
_layers:
10054: e
options: {icon: e, title: "", alt: "", clickable: true, draggable: false, …}
_latlng: o.LatLng {lat: 41.5437192894942, lng: 2.44910020452627}
_initHooksCalled: true
properties: undefined
feature: {properties: {…}}
_leaflet_id: 10054
_leaflet_events: {click_idx: {…}, click_idx_len: 1, dblclick_idx: {…}, dblclick_idx_len: 1, mouseover_idx: {…}, …}
__proto__: o.Class
10056: e
options: {icon: e, title: "", alt: "", clickable: true, draggable: false, …}
_latlng: o.LatLng {lat: 41.541313548773, lng: 2.44995964811972}
_initHooksCalled: true
properties: undefined
feature: {properties: {…}}
_leaflet_id: 10056
_leaflet_events: {click_idx: {…}, click_idx_len: 1, dblclick_idx: {…}, dblclick_idx_len: 1, mouseover_idx: {…}, …}
__proto__: o.Class`

Methods I use to add and remove layers and clusters:

function addAreaLayer(prop) {
        console.log("AFEGINT CLUSTERS "+ prop);

        active_area_list[prop]=baseLayers[prop];
        markers.addLayer(baseLayers[prop]);

        map.addLayer(markers);
        refreshActive();

};

function removeAreaLayer(prop) {
    if(map.hasLayer(baseLayers[prop])||active_area_list.hasOwnProperty(prop)){
        console.log("ELIMINANT CLUSTERS");

        delete active_area_list[prop];
        markers.removeLayer(baseLayers[prop]);

        map.removeLayer(markers);
        refreshActive();
    }
};

The markers doesn't show on the map, but if I activate the layer on map, the loaded clusters shows on map, but not dynamically by changing the bounds. Can you provide an example or explain what may be the solution?

Website where I'm working as an example of the problem:
 https://serveisweb.mataro.cat/visorSIG/images/RecursosVisor/SSIT/minisite/leaflet/MARC/visor/?perfil=patrimoni_arquitectonic


texts: Refresh layer with joined tables from watched files

texts:
I have a shapefile joined to data from a CSV file that is being watched by QGIS.

If the CSV file is updated, the attributes in the CSV layer are updated as expected, but the data in the joined attributes of the shapefile is not.

This seems like it should be automatic; why would the join still be pulling old data if the CSV layer is showing the updated data?   Am I missing a setting somewhere to make this update happen?

The only way I've been able to find to have the joined attributes show the update is if I close & reopen QGIS.  (I'm currently working in v3.4)


texts: ArcMap is telling our users that data from a certain feature dataset needs to be registered as versioned to edit.  ArcCatalog is showing us that the feature dataset in the SDE is not registered as versioned.  My question : How does ArcGIS know that the data should be versioned?


texts: Why is i.cluster not recognising my input imagery subgroup?

texts:
I'm trying to do an unsupervised image classification in Grass 7.8.  Here is what I have done so far:

Imported my aerial photography (GeoTIFF)
Set the region
Recreated the full colour image using r.composite
Recombined the image layers using i.group (created a subgroup with just the r, g, b bands)

Here is my problem: when I now try the unsupervised classification process with i.cluster I select the subgroup image that I just created  using i.group. But when I click "Run" I get the following Error message "Error in i.cluster Parameter 'subgroup'(name of input imagery subgroup) is missing".
Does anyone know what I'm doing wrong?


texts:
I have this code:

def backSDE(ActiveSDE):
    arcpy.env.overwriteOutput = True
    arcpy.env.workspace = ActiveSDE
    fcs = arcpy.ListFeatureClasses()
    print fcs

#Databases
ActiveSDE = r"C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\SERVER_SQLEXPRESS.gds\Client_Survey"
backSDE(ActiveSDE)

The database only has feature classes inside of it.  No datasets.  When I run this I get "None" in the python interpreter.  There should be 10 feature classes.  What am I doing wrong?  I'm running ArcMap 10.1

UPDATED (2)

I now have this code

def listFcsInGDB():
    ''' set your arcpy.env.workspace to a gdb before calling '''
    for fds in arcpy.ListDatasets('','feature') + ['']:
        for fc in arcpy.ListFeatureClasses('','',fds):
            yield os.path.join(arcpy.env.workspace, fds, fc)

#Set workspace environment to geodatabase
arcpy.env.workspace = r"C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\PETROPAM_SQLEXPRESS.gds\Client_Surve"

#Get list of feature classes in geodatabase
#FCs = arcpy.ListFeatureClasses()
FCs = listFcsInGDB()

 #Loop through feature classes in list
for FC in FCs:

     #Print the feature class name
     arcpy.AddMessage(FC)

This is the error message it gives me on this line from above:
    for fds in arcpy.ListDatasets('','feature') + ['']:

UPDATED (3):

Here is my updated code that lists the feature classes from the database connection file.

import arcpy
import os

arcpy.env.workspace = r"C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\PFS_Survey.sde"
def listFcsInGDB():
    ''' set your arcpy.env.workspace to a gdb before calling '''
    for fds in arcpy.ListDatasets('','feature') + ['']:
        for fc in arcpy.ListFeatureClasses('','',fds):
            yield os.path.join(arcpy.env.workspace, fds, fc)

#Set workspace environment to geodatabase
arcpy.env.workspace = r"C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\PFS_Survey.sde"

#Get list of feature classes in geodatabase
#FCs = arcpy.ListFeatureClasses()
FCs = listFcsInGDB()

 #Loop through feature classes in list
for FC in FCs:

     #Print the feature class name
     arcpy.AddMessage(FC)
     print FC
     arcpy.CopyFeatures_management(FC, "C:/Users/me/Desktop/" + os.path.basename(FC).split('.')[0] + ".shp")

It outputs this:

C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\Survey.sde\Client_Survey.DBO.CO_Arapahoe
C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\Survey.sde\Client_Survey.DBO.CO_Montezuma
C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\Survey.sde\Client_Survey.DBO.CO_Weld
C:\Users\me\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\Survey.sde\Client_Survey.DBO.Wy_Converse

The problem now is trying to copy it out of the list.  I've updated my code to how I want to copy it, but the ListFeature classes has two decimals in the code.

UPDATE 4

I'm starting a new question that is somewhat based off of this one, but different enough to be in a new page.  I will link once it's up.

Here's the link to the question about copying several sde connections.
Backup several SDE database connections


texts: Function or method in GEE that allows to fill sinks in Digital Elevation Models (DEM)

texts: I would like to know if there is any function or method in Google Earth Engine (GEE) that allows to fill sinks in Digital Elevation Models (DEM) like the algorithms available in traditional GIS softwares that use the D8 model.


texts: LiDAR (.las) to bare-earth DEM: average points spacing and cell size?

texts:
I am using ArcGIS 10.1 to generate a DEM (bare earth) from LAS dataset. To accomplish this, I use Las Dataset to Raster tool as described in Esri website.

I have read here and there in the web that the the cell size of the output raster should be set at a value corresponding to the average LAS points spacing; somewhere else, they suggest to use a value that is three times the average point spacing.

Now, since I wish to produce a bare earth DEM (i.e., with only ground returns), I am wondering if the cell size has to be set at the average point spacing of the whole LAS dataset (i.e., comprising non-ground, ground, etc), or exclusively of the ground points?

On the basis of the ArcGIS's point file information tool (example), I managed to get statistics for each LiDAR class. The average spacing (whole dataset) is 1.5 m, while for ground points only is 0.67.

Any insight into the issue?


texts: Automatically add value from another attribute column to "legend" in QGIS?

texts:
I am creating a QGIS map of events for which I have the date (year) and also rather long event names. I would like to show only the dates in the actual map and display the names in a legend. I am also using the dates to colour the points by category (see screenshot below), so "legend" shows me the dates as default values.

Instead of adding the new values manually in each row, I would like to use a function to automatically add the values from my "name" attribute column. Is that possible? I could not find what I was looking for in the QGIS-3 documentation, so I am grateful for your advice.


texts: Get intersection points between multipolygon and multiline

texts:
I am trying to get the intersection points between the boundary of a city (stored in the database as a MultiPolygon) and the OSM roads (stored as Multiline string).

I am using ST_intersect(city.geom, roads.geom) as a where clause but I am not getting any results.

If I open the tables in QGis I see that there are intersections. Do I need to transform any of the geom ?


texts: Trying to create a Cartogram and the st_transform function does not work

texts:
I'm following a tutorial to create a cartogram and the line to do so creates an error, as does the solution to fix the error which is to transform the data using st_transform(). I've looked online for similar issues/fixes and can't seem to find much.

Initial Code:

library(cartogram)
library(maptools)     # world boundaries coordinates

data(wrld_simpl)
afr=wrld_simpl[wrld_simpl$REGION==2,]

plot(afr)
afr@data

#Line causing problems
afr_cartogram <- cartogram_cont(afr, "POP2005", itermax=7)

Initial Error:

> afr_cartogram <- cartogram_cont(afr, "POP2005", itermax=7)
Error: Using an unprojected map. This function does not give correct centroids and distances for longitude/latitude data:
Use "st_transform()" to transform coordinates to another projection.

Code trying to fix previous error by calling st_Transform():

afr_proj <- st_transform(afr, "+init=epsg:4326")
afr_cartogram <- cartogram_cont(afr_proj, "POP2005", itermax=7)

New Error:

> afr_proj <- st_transform(afr, "+init=epsg:4326")
Error in st_transform(afr, "+init=epsg:4326") :
  could not find function "st_transform"
> afr_cartogram <- cartogram_cont(afr_proj, "POP2005", itermax=7)
Error in cartogram_cont(afr_proj, "POP2005", itermax = 7) :
  object 'afr_proj' not found


texts: Copying features to dataset in loop gives ERROR 000210 from ArcPy?

texts:
I'm trying to write script, that looks at the feature classes in a geodatabase and depending on what the feature class ends with, it creates a new dataset in the geodatabase and copies the file into the new dataset. I have 500 feature classes all ending with a value in the range '1 - 22' therefore the script should create 22 datasets and copies all files ending with the same number into its respective dataset. However, I keep getting an ERROR 000210: Failed to execute (CopyFeatures).This is the code I have so far:

import arcpy
import os
arcpy.env.workspace = "C:\Users\Daimon Nurse\Desktop\DFMPROJECT\DFMPROJECT.gdb"
arcpy.env.overwriteOutput=True
features = arcpy.ListFeatureClasses()
i = 1
while i <= 22:
    out_dataset_path = "C:\Users\Daimon Nurse\Desktop\DFMPROJECT\DFMPROJECT.gdb"
    out_name = "Zone{0}".format(i)
    arcpy.CreateFeatureDataset_management(out_dataset_path, out_name, sr)
    for feature in features:
        outfile = os.path.join(out_dataset_path,out_name)
        if feature.endswith('n{0}'.format(i)):
            outfile1 = os.path.join(outfile,feature)
            arcpy.CopyFeatures_management(feature, outfile1)
    i = i + 1

del i, outfile, feature, features, out_name, out_dataset_path, inlayer


texts: Updating column in layer with NULL currently to nearby value in QGIS?

texts:
Currently I have a layer with construction year attributes in it. Not all of these are filled out and I figure it is reasonable to update them with the construction year of the nearby elements as they would all be constructed around the same time typically.

How do I do this given that I would be taking values from the same table and assigning them using intersecting/crossing/contains?

I am using QGIS 2.6.0 Brighton.


texts: Mapping composite over many years in an image collection

texts:
I am looking to make an image collection that holds the 95th percentile of NDVI for each year  from 1985 - 2011 for the Landsat 5 record. I am able to create the 95th percentile composites, but I would like to avoid copy and pasting to change the year. This is what a single year would look like. How do I easily loop/map over each year to create an image collection that holds finalmask_1985, finalmask_1986, finalmask_1987.... and so on?
//////////////////////////////////////////////////////////////////////////
//1985
//////////////////////////////////////////////////////////////////////////
//Load Landsat 5 image collection
var collection_1985= ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('1985-01-01','1985-12-31')
  .map(addQualityBand);
//Isolate NDVI and time band
var only_NDVI_1985 = collection_1985.select('nd','system:time_start');
//Create image counts mask
var count_1985 = only_NDVI_1985.select('nd').reduce(ee.Reducer.count());
var countmask_1985 = count_1985.select('nd_count').gte(10);
//Make 95th percentile NDVI composite
var ninety_five_percent_1985 = only_NDVI_1985.reduce(ee.Reducer.percentile([95]));
//Create NDVI <.1 mask
var mask95_1985 = ninety_five_percent_1985.select('nd_p95').gte(0.1);
//Apply masks
var lte95_1985 = ninety_five_percent_1985.updateMask(mask95_1985);
var finalmask_1985 = lte95_1985.updateMask(countmask_1985);

Here is the link to the full script: https://code.earthengine.google.com/3674fa123d228d3a5a17b0fbcdeb3f60


texts: Why does only one of the following queries throw an error?

texts:
There are 2 SQL queries:

SELECT wkb_geometry FROM building_polygon WHERE ST_Contains(wkb_geometry, ST_GeomFromText('POINT (29.4 60.1)'));
SELECT wkb_geometry FROM building_polygon WHERE ST_Contains(wkb_geometry, ST_GeomFromText('POINT (30.32670 59.92843)'));

The first query works fine but returns 0 rows. But the second query throws an error Operation on mixed SRID geometries. It is not problem to point SRID for second query but it is not clear for me why the first query doesn't returns the same error. Is it expected behavior?


texts: How to change the display order of proportional circles in QGIS 2.10?

texts:
I'm making a map with proportional circles based on point data. I want the biggest circles to be displayed below the others and the tiniest circles above. Unfortunately the default view is different:

Is it possible to change the display setting of proportional circles in QGIS?


texts:
I'm currently working on the implementation of a Geospatial Mapping web application.  I/We  need this to be a WMS SOAP Client,  rather than the usual WMS requests via HTTP/REST.

In order to facilitate this effort I’m thinking that we could do one of the following:

1: Is there a way ( plugin, configuration) to enable GeoServer to respond to WMS SOAP requests ?

2: Does anyone know of  a  WMS SOAP Server which is open / public / internet  accessible ?

ALL OF THE URLs BELOW  Responded With

internal error - server connection terminated
Description: internal error - server connection terminated

http://193.232.117.144/cgi-bin/mapserv?map=/var/www/html/mapFiles/maingis.map&request=GetCapabilities

http://132.156.97.59/cgi-bin/canmin_en-ca_ows?request=GetCapabilities

-

-

-


texts: What does first requirement for LISA statistic mean?

texts:
According to Anselin, a LISA statistic have to satisfy two requirements:

The LISA for each observation gives an indicator of the extent of significant spatial clustering of similar values around that observation.
The sum of LISAs for all observations is proportional to a global indicator of spatial association.

The second requirement is very straightforward but I cannot understand mathematically the first requirement.

For example, he said in his paper:

  The L_i should be such that it is possible to infer the statistical
  significance of the patter of spatial association at location i. More
  formally, this requires the operationalization of a statement such as
  Prob[L_i>δ_i]≤α_i...

Is this the first requirement?


texts: Google Earth Engine: How to avoid gaps in date-wise image mosaicks?

texts:
This question is a follow-up to the excellent answer of Rodrigo E. Principe here: Mosaicking a Image Collection by Date (day) in Google Earth Engine

I am trying to apply the suggested code to gather a time series of mosaics covering the state of Bavaria. But when I display the first element of the resulting collection, I see a big coverage gap. Here's the code: https://code.earthengine.google.com/01c7ca48b7710b7ddddaa564ecc12331

Can somebody explain why this happens? And is there a way to avoid this, ideally without destroying the temporal resolution? If not, I'd also be good with a solution that combines images from a time span to fill the gap (even though the time span should be as small as possible).


texts: With schema, function postgis_full_version() does not exist

texts:
A rails application is configured to run with postgis.  Database.yml is configured as follows:

development:
  adapter: postgis
  database: app_development
  username: simpleton
  schema_search_path: public, postgis

for the schema search path, the attempt is to ensure that postgis is not installed twice by enacting as follows
app_development=#

DROP EXTENSION PostGIS;

CREATE SCHEMA postgis;
CREATE EXTENSION PostGIS WITH SCHEMA postgis;
GRANT ALL ON postgis.geometry_columns TO PUBLIC;
GRANT ALL ON postgis.spatial_ref_sys TO PUBLIC;

which generates the error

ERROR:  function postgis_full_version() does not exist
LINE 1: SELECT PostGIS_full_version();
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Oddly, if the configuration is without a schema, i.e.

schema_search_path: public, postgis
CREATE EXTENSION PostGIS;

the SELECT command runs and the database migration can occur.

Why is this behaviour apparently selective?


texts: Mapping proportional Data for two group on single Map?

texts:
I have no experience in mapping but for my research purpose I required to map certain data on Map of certain country. Data is like this :

State     Proportion of X     Proportion of Y     Rest
A         50                  30                  20
B         79                  20                  1
C         95                  4                   1
.
.
.

Such that there are total 30 states and for each state I required to show proportion of X and Y on the same map for comparison purpose so that reader could understand the where proportion of X is high or low in contrast to Y.
One of my friend told me about software ArcMap but another told me it is not possible to plot both the proportion on the same MAP.
I have access to ArcMap or any other that would be recommended if it is open source. I have no prior knowledge in mapping in any software.

I want to know is it possible to show such data on MAP in the way what I want ? And which software is most appropriate for it [please be restricted to ArcMap or  open source software]

You may recommend me tutorial video for this purpose if providing entire solution is not possible.

I want something like this :

Here is Zoom view


texts: How to get points in bounding box by degrees for map

texts:
I have some concerns about working of BOX3D PostGIS function. I tried to find places in square -158.721365889063 66.09275164148659, 154.608712235937 52.56807548751039(degrees). But this where statement

    some_place.point && SetSRID('BOX3D(-158.721365889063 66.09275164148659,
    154.608712235937 52.56807548751039)'::box3d,4326)

returns places in the biggest of two possible squares.

When I change order of angles or use 'xmin ymin,xmax ymax' order, I got the same result.

Also I used

     ST_Contains( SetSRID( ST_MakeBox2D(
      ST_Point(143.930001298437,66.82980906973742),
      ST_Point(-145.757498701563,51.37699768868392) ), 4326 ), some_place.point)

with the same result.

How can I do it for some specific order of angle's points?

P.S.
I got bounds from google map api.


texts: Trying to link ZIP code to core-based statistical area

texts:
I have a dataset downloaded from the IRS reporting AGI by ZIP code (link: irs.gov). However, I'm interested in tying each ZIP code to statistical areas (whether it's [a] primary statistical area, [b] combined statistical area, [c] core-based statistical area, [d] metropolitan statistical area, or [e] micropolitan statistical area). Having all of them in one dataset would be ideal, but I'd settle for just CBSA.

I found a great "crosswalk" dataset on the HUDUser website (link: huduser.gov) that does exactly this--connects ZIP code to CBSA. I merged the two files in SAS on ZIP code and voila, most of my problems were solved. But I realized that there were some blaring issues with the merged dataset; for instance, Los Angeles was characterized by different CBSA codes in each file, amongst other errors. I went back to the download sites (both the AGI by ZIP code dataset from the IRS as well as the HUDUser ZIP code to CBSA dataset) and noticed that they were different dates--2014 and 2016Q3 (HUDUser datasets are updated quarterly), respectively. Of course, this could have been the issue, so I tried a date that was closer to 2014 (2013Q4, which was updated 2013-12-22). Fewer errors this time, but still many errors.

Can anyone weigh in on this? How can I ensure that there is no mis-match between CBSA coding on my IRS dataset and my HUDUser dataset? Or, is there another "long" dataset that contains all statistical areas by ZIP code, where each valid ZIP code is a row in the dataset?


texts:
I have this raw text value extracted from Postgres geometry(point) type 0101000020E6100000FAFFEFE0FB444B40672BE61C924E3840
I was wondering if there is any way to get the lat,long values from it without having to insert it into db again and using the ST_X, ST_Y functions!


texts:
In QGIS 1.8 I used the openstreetmap plugin to request data from the server, or open a .osm XML file, which produced 3 layers (polygons, polylines and points) in QGIS.

This was great (give or take the 64bit bug which highly affects the freshly mapped area I work on), but how do I do that in QGIS 2.x?

I can't find any updated documentation about it.

In the vector menu, I've tried importing from osm servers, I get a .osm file, but then I'm stuck. There's a menu option for loading an XML file, but it produces a spatialite db, which apparently only has non-geometry tables when I try to add it as a spatialite layer.
If I open the attribute table, I can see node ids, ways and so on, but I can't figure out how to actually make use of them.

I've also tried adding a vector layer from a .osm XML file. I get a gpsbabel window that shows nothing, then finally an invalid data source message in QGIS.


texts: PyQGIS adding polygon features to temporary vector layer results in empty layer

texts:
I'm trying to extract certain features from a layer and add them to a temporary layer to run an algorithm on in a later stage to limit processing time. Code below selects the correct features and adds them to the temporary layer, but when I check the layer, it shows no geometry and no attribute information even though I ask to copy these. It only gives the correct number of intersecting features (see image below). I've searched extensively for a solution for this, but can't find an appropriate one.
What am I missing?
# create temporary layer
vl = QgsVectorLayer('Polygon', 'temporary_polygon', 'memory')
pr = vl.dataProvider()

idx = joinOutput['OUTPUT'].fields().indexFromName('OPPERVL')
for perceel in joinOutput['OUTPUT'].getFeatures():
    # to run code faster, limit to a features within a certain area interval
    if perceel.attributes()[idx] > thisArea*0.9 and perceel.attributes()[idx] < thisArea*1.1:
        for gw in GWPov.getFeatures():
            if perceel.geometry().intersects(gw.geometry()):
                # add intersecting feature to temporary layer
                feat = QgsFeature()
                feat.setGeometry(perceel.geometry())
                feat.setAttributes(perceel.attributes())
                pr.addFeature(feat)
                print('Test')

            else:
                continue
    else:
        continue

# control check, showing empty geometry and attribute table
QgsProject.instance().addMapLayers([vl])


texts: Adding background image to ArcGIS online Popup

texts:
I am trying to make a custom popup for a storymap template by changing the CSS.  The background colors can be changed but is there a way to add an image to the main area with the attributes instead of a color background.   I have tried

.esriPopup.chrome .titlePane,
.dj_ie7 .esriPopup.chrome .titlePane .title{
    back-ground image: URL (/images/image.png)
}

and

 .esriPopup.chrome .titlePane {
    back-ground image: URL (/images/image.png)
}

Nothing has worked.  Is this possible and if so where would I put the back - ground image code?


texts: Compare features in PostGIS with updating attributes

texts:
I have two polygonal layers: "existing" and "modified".

Attribute table of "existing" and "modified" layers have same column "feature_id" - number of feature. Please help me create query which write attributes 't' (if features matches) and 'f' (if features not matches) in 'compare_result' column of "modified" layer using ST_Equals function.
Now I receive attributes 't' (if features matches) and 'f' (if features not matches) in pgAdmin result window using query:

SELECT  ST_Equals(modified.the_geom, existing.the_geom) AS compare_result FROM  modified, existing
WHERE modified.feature_id=existing.feature_id

but I don't can create query? which write result in attribute table. Thanks!

EDIT:
Solution without creating additional layer:

UPDATE
    modified
SET
    compare_result = compare_result_temp

FROM
    (
        SELECT
r."feature_id" AS idmodified,
ST_Equals(r.the_geom, m.the_geom) AS compare_result_temp
        FROM
     "modified" AS r,
     "existing" AS m
WHERE r."feature_id" = m."feature_id"
    ) AS foo
 WHERE
 foo.idmodified = modified."feature_id";


texts: QGIS crashed when click on anything in plugin manager

texts:
When I click on anything inside of the plugin manager QGIS crashes immediately. (QGIS 3.0.3, arch linux)
Running from the command line this error shows up on crashing:

free(): invalid pointer
[1]    29800 abort (core dumped)  qgis


texts: I am trying to simply get the value of each cell in my raster using the Raster To ASCII tool in ArcGIS Desktop.  I have absolutely no idea what I am doing wrong. The output I get is all -9999 when the input rasters are float with a range of values of about -3 to 35. If I use Raster To Point then I get the actual values to the points, but I need these data in matrix/table form.  What the heck am I doing wrong, I have tried everything that I can possibly think of!


texts: Merge polyline feature at vertex intersection based on segment count

texts:
I have a Polyline layer in QGIS which is split into multiple segments using the explode lines tool (the black lines in the image). There are nodes created from the 'extract vertices' tool (the red circles in the image). I want to merge the segments but keep 3-way or greater intersections split. Dead ends can be ignored but must be kept.
When splitting the features I could not find a way to split exclusively at vertex intersections of 3 or more. I then decided to explode the line and try merging intersections with less than 3 segments, but I'm stuck again. Ideally this method can be automated.


texts: Making digitization more detailed afterwards with QGIS?

texts:
I digitized already two areas, bordering each other on one side. Then I realized that more detail was needed!

What do you suggest as a fast way to do so?

Adding each time new dots/vertices takes so much time.
But, adding a new polygon and merge both (as I used to do for a not bordering polygon) does not seem that easy in this case.
Using the split feature causes an invalid area which I cannot merge with the other area afterwards.


texts:
I'm trying to find any documentation to help me understand and use the leaflet-geopackage plugin.
what I want to do is use the full potential of the .gpkg file as a database and query only part of the layer.
is that even possible? or do I have to just hide the parts I don't want to show?

Clarification, I'm only talking about feature layers, not tiles.
loaded with the example code from their github:

var rivers = L.geoPackageFeatureLayer([], {
     geoPackageUrl: 'http://ngageoint.github.io/GeoPackage/examples/rivers.gpkg',
     layerName: 'rivers',
     style: {color: 'green'}
    }).addTo(map);


texts: Excluding point from Nearest neighbor search once it's been matched using R?

texts:
I have two lists of addresses, List 1 and List 2. I want to find the nearest neighbors of the list--but without assigning more than one point to any point. So say point 1A is the nearest neighbor of points 2A and 2B, but point 2A is closer. Then I want point 1A to match with point 2A, and point 2B to match with its own next nearest neighbor. Basically, every point is paired and there are no 3+ point clusters.

Is this possible to do using R?


texts: Losing mouse focus on map in safari browser

texts:
I am facing problem when I use SAFARI browser for GIS WEB Application. I have select tool in my application to select the features on the map. When i use Safari the point focus is not correct it drags to some other location.

I have tried the same functionality in other browsers like IE, Chrome and Mozilla Firefox. It works perfectly.

Any suggestion what additional needs to be done to make it work in Safari Browser?

Also is Safari recommended for GIS Web Application?

The snapshot1 is of Safari web browser where the mouse focus is lost and it goes out.

Whereas the actual result in other browser is in snapshot2.


texts:
I have a .shp-file of some buildings with a sheet in an attribute table of their height. I want to extrude them using a Qgis2threejs plugin, but since a sheet of height is shown as a type of QString, the plugin can't recognize it and the buildings can't be extruded.

I have tried to export a database as a .csv-file, but the Excel is locked and it can't be modified.

Any alternative approaches?


texts: Display Network analysis route service in India

texts:
I am working on ArcGIS JavaScript API, I need to display route from one point to another.
Thats why I need any Network Route Service for India. Is there any?


texts: Google Maps API for iPhone and Android application

texts:
I need to develop map with Google Maps API with markers etc. That map would need to be integrated into iPhone and Android application. That integration will be done by another person.

My question is: is it possible to develop map, define all coding with Google Maps API and then just give that coding to the person which developed application to integrate it into iPhone and Android application? I don't know how developing mobile applications works so I don't know if that's ok.

Thank you.


texts: Convert Polygon exterior ring to points and linestring to points in SQL

texts:
I am try create materialized views in Oracle listing a polygon and linestring vertices with their longitude and latitude in separate columns. So I know I would need to count the vertices in each polygon or linesegment but I am not having any luck on either task.

I am trying to get a loop to work to count the vertices in each object but it's not working.

for nn in 1..(nn-1) loop
select OBJECTID,
SDE.st_x (SDE.st_pointn (sde.st_exteriorring (SHAPE), nn)) as X,
SDE.st_y (SDE.st_pointn (sde.st_exteriorring (SHAPE), nn)) as Y
from
(SELECT sde.st_numpoints (shape) as nn from dual, SDW.trpBridge_test), SDW.trpBridge_test
order by objectid;
end loop;


texts: Get feature from ol.interaction.Draw on drawend

texts:
My problem seems simple but I cannot find any questions that successfully solve it. I want to draw a circle using ol.interaction.Draw, and access the feature that is created immediately.

My code as of now is:

var radiusSelectInteraction = new ol.interaction.Draw({
    source: circleLayer.getSource(),
    type: 'Circle',
    features: circleFeatures
});

radiusSelectInteraction.on('drawend', function () {
    circleFeatures.forEach(function(feature) {
        ....
    });
});

My problem with this is that it doesn't give me the feature that is currently being drawn. If I don't clear the vector, it will show the first circle on the drawend event of the second one I draw. I only plan on allowing 1 circle at a time so I need to access the feature on end of drawing or on end of feature add.

I tried something like this instead of drawend:

radiusSelectInteraction.on('drawstart', function () {
    circleLayer.getSource().once('addFeature', handleAddFeature);
});

function handleAddFeature(evt) {
    var feature = evt.feature;
    ....
}

but it was never triggering the handleAddFeature event like it was supposed to. Did I do it wrong or is there a better way to achieve this?


texts:
At which (xth) decimal place of the coordinate values QGIS defines any set of points is duplicate?

I had been thinking QGIS handles ~15th decimal place; but it was the limitation only because I had been working on Shapefiles mainly.

However, the limit does not apply when I choose another data source, like a temporally scratch layer.

Using a dummy data like below, I made a quick test by two tools, to find the smallest value they can identify differences of coordinates:

QGIS geoprocessing: Delete duplicate geometries
MMQGIS plugin: Delete Duplicate Geometries

Output by QGIS Delete duplicate geometries is the same as the above input data (all 20 records were preserved), so it may mean QGIS thinks they are all distinct. Does this limit go beyond 1e-29 (or 1e-30) seen in this small test?

Just as comparison, MMQGIS Delete Duplicate Geometries produced below. It seems MMQGIS sets the 16th or 17th decimal places as the limit.

[Edit]

I am afraid I had not been clear what was the central part of my question. My aim to understand the relationship between coordinate values and the duplicates/overlapping, that is as stated in the first paragraph of my question. Hopefully this kind of knowledge can help us easily control overlapping features by tweaking field calculator expressions.

However, underlying issue I had been trying to focus was that the decimal places QGIS recognizes points (nodes) as overlapping seem to vary, depending on the data source.

If we use Shapefiles for our layer, QGIS handles ~15th decimal place and smaller difference (at 16th or 17th) is not recognized... MMQGIS is also affected by this limitation (to my eyes). If we create a temporally scratch layer or DB layers this limitation goes well beyond 30th-ish? This change in behavior makes me wonder.


texts:
I know the question has been asked before, but I get an error trying to clip my nc file with a shapefile with the method used in this post.
My error is DimensionError: x dimension not found. 'set_spatial_dims()' can address this.
My code is
path = 'a/given/path/'
import geopandas
import rioxarray
import xarray
from shapely.geometry import mapping

nc = xarray.open_dataarray(path + 'FR.nc')
nc.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
nc.rio.write_crs("epsg:4326", inplace=True)
shp = geopandas.read_file('AK.shp', crs="epsg:4326")

clipped = chen.rio.clip(AK_Shape.geometry.apply(mapping), AK_Shape.crs, drop=False)

The nc file and the shape file has the following dimension and information
nc: <xarray.DataArray 'FR' (lat: 180, lon: 360, model: 15, regr: 5)>
shp: Index(['AREA', 'PERIMETER', 'STATESP010', 'STATE', 'STATE_FIPS', 'ORDER_ADM',
       'MONTH_ADM', 'DAY_ADM', 'YEAR_ADM', 'geometry'],
      dtype='object')

Can I do a work around and then get the cut-out version of my nc file with the shapefile?


texts: How to access map-box fill-extrusion layers to apply hover effect/popups

texts:
I am looking for a way to create a hover effect on a fill-extrusion layers

I have replicated the demo at: https://www.mapbox.com/mapbox-gl-js/example/hover-styles/   [this is on flat shapes]

The I have also built a demo using fill-extrusion layers layers where the height and color is set based on an input set

Both the above samples works fine

In my fill-extrusion project, I cant seem to crack the code on hovering over the layer - when testing with the simple function:

map.on("mouseenter", "state-fills", function() {
    // Change the cursor style as a UI indicator.
    map.getCanvas().style.cursor = 'pointer';
});

I don't see the the cursor changing.

This leads me to believe that the "mouseenter" on layer with id: state-fills is not the correct way to achieve this?

Ultimately I would like to create a hover effect with a popup, but for now, If i can just figure out why my "mouseenter" doesn't do anything, Im sure I will be going in the right direction

I have also tried replacing mouseenter with mousemove - but to no avail

Here is the more complete map init:

map.on('load', function () {
  // Include the Province Shapefile
  map.addSource('mySrcName', {
    type: 'vector',
    url: 'mapbox://pantoeflou.82pwfdao'
  });

  //Create a layer from the source data -- Works Fine
  map.addLayer({
    source: 'mySrcName',
    'source-layer': '3_Fields-83vr21',
    id: 'state-fills',
    type: 'fill-extrusion',
    paint: {
      'fill-extrusion-opacity': 0.75,
      'fill-extrusion-color': {
        property: 'PROV_ID',
        type: 'categorical',
        stops: colorStops('Val') // Data attbt
      },
      'fill-extrusion-height': {
        property: 'PROV_ID',
        type: 'categorical',
        stops: hightStops('Val2') // Data attbt
      }
    }
  });

    // !! Does Nothing !!
    // When the user moves their mouse over the states-fill layer, we'll update the filter in
    // the state-fills-hover layer to only show the matching state, thus making a hover effect.
    map.on("mouseenter", "state-fills", function() {
        // Change the cursor style as a UI indicator.
        map.getCanvas().style.cursor = 'pointer';
    });

});


texts: SyntaxError: non-keyword arg after keyword arg python

texts:
I'am working a project about hot clouds volcano mountain simulation with Titan2D v4.0.0 but I got error when I try to run the simulation.py
sim.setGIS(
   gis_format='GIS_GRASS',
   gis_main='/home/dinar/Downloads/titan2d-v4.0.0/grass',
   gis_sub='newLocation',
   gis_mapset='PERMANENT',
   gis_map='merapi8m.grass',
   gis_vector=None,
   region_limits==(417661e+05, 9143861e+06, 4611590e+05, 9176061e+06)
)

it said non-keyword arg after keyword arg python


texts: Distance between points along polyline ArcGIS for Desktop?

texts:
I'm using ArcGIS 10.3 (Desktop), and I'm trying to measure the distance between points following a path determined by a set of polylines.  Hand measurement is not a viable option-- there are ~7500 points, and my polyline set is the global coastline.  The ultimate goal is to get the spacing of these points along the polylines, so that I can analyze this with respect to a few other variables.

I can read python-- with a bit of effort-- but not write or edit it, so code-based solutions will not be particularly useful to me.

I have tried using the linear referencing tools-- I snapped the points to the lines, then transformed the polylines to routes and used 'locate features'-- but because my points are not (entirely) ordered along the path, using the difference between  linear reference locations of subsequent points yields some fairly unlikely answers.
In addition, spot checks of (hand-measured, using the built-in path length measuring tool) point-to-point along-polyline distance yield different answers than given by the corresponding linear references for the point pair. (Some of the answers were close enough that they could very well be measurement error on my part, but others were nearly 20%, despite my having zoomed in and snapped the measurement cursor to the polyline, proceeding vertex-to-vertex.)

Sorting first by linear distance seems to fix the first problem, but in light of the second one, I remain skeptical of this process.

I then tried simply splitting the polylines into segments ('split line at point') using the (snapped) points as the dividers, with the eventual goal of finding the segment length with the calculate geometry feature.  This produced gibberish. I'm not sure what criteria ArcGIS used to split the lines, but it wasn't the point set I provided.  Segments in the resulting shapefile ran right through tens of points, only to terminate at a place with no point. (Also, there were segments utterly empty of points...you get the gist.) I reiterate that the pointset was specifically snapped to the polylines, so this problem shouldn't be caused by excessive distance between the point and the polyline.

I'm at a bit of a loss, here: could someone point out where I might have gone astray and/or suggest a solution?  Or, failing that, a new method that I could try?  It seems like a fairly basic (not necessarily simple, but fundamental) analysis to do with geographic data, and I don't understand why it has presented so many difficulties.


texts: How can I join a buffer layer to the layer it was created from?

texts:
I have two relevant layers for this problem:

Real estate parcels (point data)
Multiple ring buffer layer - contains multiple buffer rings around each particular real estate parcel. (0.5, 1, 1.5, and 2 miles).

There are 2900 real estate parcels, each one a multiple ring buffer layer containing 4 quartiles. The buffer layer table has 2900*4 entries since each parcel has 4 buffer layers.

I need to join the buffer layers with the parcel they were created from, and there doesn't seem to be an easy way. Any tips on a useful tool?


texts: Working with remote QGIS server project files

texts:
I went into some troubles when making my first steps with QGIS server and I could not find an easy way to go.
Some context
My projects files (qgs/qgz) are on a remote server (Ubuntu based, headless, i.e. no GUI), where they are exposed to the web, mainly for serving raster files through the OGC WM(T)S protocol and some vector data from PostGIS to many depending applications in a larger infrastructure.
Whenever I want to change/update something in a project, I have to copy the project file over SSH to my local computer, work on it, then push it back again on the server.
This is not optimal, mainly for 3 reasons;

Version mismatch: my local version of QGIS is more recent than on the server. Hence, I may get into troubles...
File paths: paths of, for example, raster data are obviously not the same on my machine than on the server. I have to tweak in either way or the other, and it's especially painful to work with extra large raster. I don't want to transfer these files locally. And trying to simulate their existence locally by small subsets is OK but it involves some manual processing and, therefore, is prone to errors.
Remote database: I've to manage to connect to the remote PostGIS using an SSH tunnel or to duplicate the database locally, which in case of a large database is not possible.

Some pseudo-solutions
With these findings, I see two possibilities (but I probably missed the obvious ones), with pros and cons for each:

Editing the raw XML: I tried to directly edit the XML based project file on the server, but it's definitely not the way to go as it becomes a really huge file; you quickly miss something or even break it. And this of course doesn't work with qgz files.

Installing a GUI: In order to spare resources, I also want to avoid the installation of a graphical interface on the server (which would naturally make it possible to remotely work on the project solving all the previously described problems, but which need a good VPN connection).

I probably don't work in the right manner (is there just one?) with remote QGIS project files which are used/served by QGIS server.
Question
So, how do you manage to work cleanly and rigorously with remote QGIS project files, are there currently recognized best practices?


texts: How to place a label below a symbol in QGIS 3.18.3

texts:
I have an airport symbol and below that I want to show the label for the airport identifier.
Because I want to keep a proportional distance between the label and the symbol at different zoom levels, I used and offset as 'Meters at scale', set the symbol to be 'offset from point', on the lower quadrant, and offset 8000 on the y axis.
With zero offset on the y, the label shows over the airport symbol. With 8000 on the y it moves up. However with -8000 the label remains in the same place. It doesn't move down as I would expect.
How can I properly make the label show below the SVG symbol, and remain there at different zoom levels?


texts: Point pattern fill with only full points in QGIS

texts:
I would like to fill a fairly complex polygon with a point pattern fill, but the border of the polygon "cuts" the circles and the result does not look as clean as I would like.

I know it is this is how it should be, but is there a way to use only full circles/points/dots? For example, showing a black dot only if its centre is contained within the polygon (and even if part of it is outside the polygon).


texts: Error adding store through REST API to GeoServer

texts:
I have problem with adding raster layer to geoserver automatically. I need, that geoserver add layer without person job, only script action and i saw that it may be done with cURL REST API but, i can't do it.
In general, I have some GeoTIFF files and i want to published them at geoserver through the terminal. I do how wrote DaveParillo. I can add workspace,

 curl -u admin:geoserver -v -XPOST -H 'Content-type: text/xml' -d '<workspaces><name>geotiff</name></workspaces>' http://192.168.255.197:8080/geoserver/rest/workspaces

answer

* Hostname was NOT found in DNS cache
*   Trying 192.168.255.197...
* Connected to 192.168.255.197 (192.168.255.197) port 8080 (#0)
* Server auth using Basic with user 'admin'
> POST /geoserver/rest/workspaces HTTP/1.1
> Authorization: Basic YWRtaW46Z2Vvc2VydmVy
> User-Agent: curl/7.35.0
> Host: 192.168.255.197:8080
> Accept: */*
> Content-type: text/xml
> Content-Length: 43
>
* upload completely sent off: 43 out of 43 bytes
< HTTP/1.1 201 Created
< Date: Thu, 10 Dec 2015 14:15:53 GMT
< Location: http://192.168.255.197:8080/geoserver/rest/workspaces/geotiff
* Server Noelios-Restlet-Engine/1.0..8 is not blacklisted
< Server: Noelios-Restlet-Engine/1.0..8
< Transfer-Encoding: chunked
<
* Connection #0 to host 192.168.255.197 left intact

but when i try add store i touch this error:

curl -u admin:geoserver -v -XPOST -H 'Content-type: text/xml' -d '<coverageStore><name>fff</name><enabled>true</enabled><type>GeoTIFF</type><url>/opt/tomcat/webapps/geoserver/data/avhrr_locale1e2.tif</url></coverageStore>' "http://192.168.255.197:8080/geoserver/rest/workspaces/geotiff/coveragestores?configure=all"

answer:

* Hostname was NOT found in DNS cache
*   Trying 192.168.255.197...
* Connected to 192.168.255.197 (192.168.255.197) port 8080 (#0)
* Server auth using Basic with user 'admin'
> POST /geoserver/rest/workspaces/geotiff/coveragestores?configure=all     HTTP/1.1
> Authorization: Basic YWRtaW46Z2Vvc2VydmVy
> User-Agent: curl/7.35.0
> Host: 192.168.255.197:8080
> Accept: */*
> Content-type: text/xml
> Content-Length: 155
>
* upload completely sent off: 155 out of 155 bytes
< HTTP/1.1 500 Internal Server Error
* Server Apache-Coyote/1.1 is not blacklisted
< Server: Apache-Coyote/1.1
< Content-Type: text/plain
< Transfer-Encoding: chunked
< Date: Thu, 10 Dec 2015 14:16:14 GMT
< Connection: close
<
* Closing connection 0
:Store must be part of a workspace

Geoserver log same error:

2015-12-10 17:15:53,801 INFO [catalog.rest] - POST workspace geotiff
2015-12-10 17:16:14,687 INFO [org.geoserver] - Loaded store 'fff',    enabled
2015-12-10 17:16:14,688 ERROR [geoserver.rest] - :Store must be part of a   workspace
2015-12-10 17:16:14,688 ERROR [geoserver.rest] -
org.geoserver.rest.RestletException: java.lang.IllegalArgumentException:   Store must be part of a workspace

What's wrong?

Geoserver v.2.8.1.


texts: Classifying a Raster with Unique Values from Existing Layer File and Dynamically Change Labels - How to Create a UniqueValues Raster Layer File?

texts:
I'm creating a python tool that needs to have the output raster added to the current map and symbolized using unique values. I have been able to accomplish this by creating a template layer file and applying it to my output raster using arcpy.mapping.UpdateLayer . However, I also need the labels for the classification scheme to dynamically update based on user inputs to the tool. I understand that this could be accomplished using the Unique Values Symbology Class functions. However, this can only be done if the symbology type of the layer is "Unique_Values".

My issue is that I can't seem to create a raster layer file that is a "Unique_Values" type. Every time I try to create a layer file with raster - unique values symbology the SymbologyType returns as "Other". This "Other" type isn't supported in the arcpy.mapping.UniqueValuesSymbology functions. Therefore, I am unable to change any of the symbology labels.

Does anyone know if there is a way to create a unique values layer file for a raster layer that has a SymbologyType of "Unique_Values"? I have a feeling that this type isn't supported, but I would love to be proved wrong!


texts: Loading Tab Reader Parameter into ArcMap ModelBuilder

texts:
In our office, we have a large set of MapInfo TAB files, and are planning to migrate those data to ArcGIS Geodatabase with the same symbology available in the TAB file.

On a individual TAB file to GDB, we can do it like:

Import data from TAB A to GDB A using Data Interoperability, by its corresponding CRS.
Then adding GDB A's featureclass to ArcMap and "load symbology from TAB file" using TAB Reader though unique data
Then save applied symbology to ".lyr" file or converting them to "Representation".

The following figure explain the workflow, not the ModelBuilder.

To convert hundreds of TAB files is taking a lot of time. So, I thought of doing with coding or ModelBuilder. Since I am very poor with Python or any other coding, I am going with ModelBuilder. The problem arises when I tried to load symbology from TAB file.

Could any one help me out in this regard?


texts: How to render mixed geometry from one table?

texts: I have a table in postgresql with mixed geometry, but it is not accessible in QGIS. I am running QGIS version 2.8.2. Is there any plugin to install for QGIS to render multiple geometry?


texts: Download opentopo maps layer in QGIS 3.2.2 and adding GPS coordinates

texts:
I would like to create a topographic map around Kinanira (Uganda) with QGIS, in which I can load my GPS coordinates from my soil samples.

In order to get a better overview, it would be useful to color the reports according to different pH values.


texts: Given an array of OSM point IDs, how can I create a path of connected roads?

texts:
I'm writing a shortest route script in Python using geospatial data from a Postgres DB with PostGIS extension and pgRouting to make the points from OSM data.

I have a path of OSM IDs by points. I got this from downloading OSM XML data from a particular region, importing the data into a database using osm2pgrouting, calculating the shortest route in my script, and using the OSM ID as the main point of identification. This all being stored in a PostGIS database.

If I had an array of, in order from 0-n for the path, the OSM IDs of all the points for the shortest route, how could I create a line to line up with the roads in OSM for the path?

I've looked at ST_MakeLine, but my worry with that is it will output like the picture below, and not connecting the actual roads. I want the path to actually map with the actual roads and not just create a line between points where a road doesn't exist.

What should I do?


texts: Start-up error message QGIS 1.8 Couldn't load PyQt4

texts:
After successfully building QGIS 1.8 "Lisboa" on Windows 7 in 32 bit mode by following the instructions mentioned here I am able to start QGIS, however I get the following startup error message:

  Couldn't load PyQt4. Python support will be disabled.

  Traceback (most recent call last):   File "", line 1, in  ImportError:
  DLL load failed: The specified procedure could not be found.

  Python version:
  2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

  QGIS version:
  1.8.0-Lisboa 'Lisboa', exported

After this message the Python console and other python related plugins are not loaded.

I have installed the latest PyQt installer of PyQt4.9.5 from riverbanks

My system details

OS: windows 7 64 bit

Build tools:

MS Visual Studio 2010
C-Make 2.8.8
Python2.7
PyQt4.9.5 with SIP
OSGeo4W with GRASS and all dependency libs

further investigation lead me to believe that there might be a mismatch of PyQt 4.9.5 that I installed and the one downloaded by OSGeo4W (PyQt4.8). However, I do not know how to make the cmake look for pyQt4.9.5 instead of PyQt4.8

I would highly appreciate any pointers and detailed help in solving this issue.

With regards

Sukesh


texts: What are the advantages / disadvantages of reporting location in decimal degrees versus degree min sec?

texts:
Many scientific articles report site location in units of degree/min/sec, but for analysis in with non-GIS software, it is easier to use decimal degrees. Is there a reason for degrees / min / sec to be the more frequently used convention? What are the advantages and disadvantages to using one or the other?

clarification While I am generally interested in the advantages and disadvantages to using one or the other, my particular interest is in the context of reporting location in a scientific paper.


texts: Walking-Driving area with Leaflet,OpenLayers

texts:
Could be possible to get a walking or driving distance area using pgr_drivindistance,pgr_alphashape and/or pgr_pointsaspolygon when an user click any location  on map with Leaflet or OpenLayers?.

I am new programming in sql and trying to get the correct result I found the first problem on how to pass source node which is integer to double ( x , y).

The second one is how to join pgr_drivindistance and pgr_alphashape so that through Geoserver get the walking-driving polygon when the user clicks(x,y) on the map.

Would have another chance to achieve that objective with another functions?.


texts: Labeling in QGIS: what does the priority bar under "position"

texts:
when you label a layer in QGIS, under the tab "position" at the very bottom, you can set the priority from low to high.

I can't find in the documentation what it does. What is beeing prioritized?

Sorry if it is in the docs, but I don't find it.


texts: Get OSM building heights using height or building:levels tags

texts:
I want to parse OSM data and to get all buildings with height or building:levels data.
I changed osmconf.ini file to read height and building:levels tags as separate columns for multipolygon layer but after all this i get an error ERROR 1: "height" not recognised as an available field. ERROR 1: SetAttributeFilter(building is not null and (height is not null or building_levels is not null)) on layer 'multipolygons' failed.
This is ogr2ogr command I used ogr2ogr -f "GPKG" -where "building is not null and (height is not null or building_levels is not null)" building_height.gpkg us-latest.osm.pbf multipolygons -progress --config OSM_CONFIG_FILE "C:\\Program Files\\QGIS 3.10\\share\\gdal\\osmconf.ini"
I have no idea why it is not working.
Here is osmconf.ini snippet for mulitpolygons:
[multipolygons]
# common attributes
# note: for multipolygons, osm_id=yes instantiates a osm_id field for the id of relations
# and a osm_way_id field for the id of closed ways. Both fields are exclusively set.
osm_id=yes
osm_version=no
osm_timestamp=no
osm_uid=no
osm_user=no
osm_changeset=no

# keys to report as OGR fields
attributes=name,type,aeroway,amenity,admin_level,barrier,boundary,building,craft,geological,historic,land_area,landuse,leisure,man_made,military,natural,office,place,shop,sport,tourism,industrial,addr:city,addr:state,addr:street,addr:postcode,addr:housenumber,addr:unit, height, building:levels
# keys that should NOT be reported in the "other_tags" field
ignore=area,created_by,converted_by,source,time,ele,note,openGeoDB:,fixme,FIXME
# uncomment to avoid creation of "other_tags" field
#other_tags=no
# uncomment to create "all_tags" field. "all_tags" and "other_tags" are exclusive
#all_tags=yes

[multilinestrings]
# common attributes
osm_id=yes
osm_version=no
osm_timestamp=no
osm_uid=no
osm_user=no
osm_changeset=no

# keys to report as OGR fields
attributes=name,type
# keys that should NOT be reported in the "other_tags" field
ignore=area,created_by,converted_by,source,time,ele,note,openGeoDB:,fixme,FIXME
# uncomment to avoid creation of "other_tags" field
#other_tags=no
# uncomment to create "all_tags" field. "all_tags" and "other_tags" are exclusive
#all_tags=yes


texts: Seeking ST_Subdivide implementation in Python?

texts:
I'd like to utilize PostGIS' ST_Subdivide operation via Python. I would like to avoid interfacing directly with PostGIS, have a database up, etc. Instead, I would like to pass a geometry object and have the function evaluated in memory.

I can see that ST_Subdivide wraps over the function lwgeom_subdivide which itself wraps a recursive method defined here.

I am looking for either:

An implementation of this algorithm in Python
An isolated version of this method that has Python bindings... or some resources that might help me go about doing this
A pattern for potentially accessing this function from Python (likely not a reasonable request)


texts: Removing layer and exporting mapbook into PDF?

texts:
I am new Python.

I need to write a script that turns off a layer (and removes it from the legend, which it would probably do automatically since that box is checked in legend properties) and then I need to export the map book.

I need to do this many times, so it'd be nice to have a script.

Heres what I have. It runs and gives me no errors and also doesn't give me a pdf...

import arcpy
from arcpy import env

'# Local Variables

Goldendale = r"\\WSBS\Company\PROJECTS\2013\12 - December\131226 - Goldendale SMP\GIS\Data\Shapefiles\goldendale_UGA.shp"

'# Remove Goldendale

mxd = arcpy.mapping.MapDocument(
    r"\\ws01\Company\PROJECTS\2012\12 - December\121201 - Klickitat SMP\GIS\MXD\InventoryMapfolio\10_KlickitatSMP_FloodplainFloodways.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
layers = arcpy.mapping.ListLayers(mxd, "*", df)
names = [Goldendale]

for layers in layers:
    if layers.name in names:
        layers.visible = False

arcpy.RefreshTOC()
arcpy.RefreshActiveView()

'# Save

mxd.save()

'# Export to Mapbook to PDF

for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
    mxd.dataDrivenPages.currentPageID = pageNum
    arcpy.mapping.ExportToPDF(mxd,r"\\ws01\Company\PROJECTS\2012\12 - December\121201 - Klickitat SMP\GIS\PDF\InventoryMapfolio\Revised_20160812\10_KlickitatSMP_FloodplainFloodways.pdf", "ALL")

What would I need to change in order for me to get a combined PDF without the layer?


texts: Problem with the correct calculate of the bisector of the angle

texts:
I have a problem with the correct calculate of the bisector of the angle. I use WGS84 projection. I have a bisector's  bearing (azimuth), coordinates of start point and line length as input. Some bisectors are determined correctly, but I have a problem with this case:

input:

start_point=[19.299301570524328, 50.35459214539403]

bearing=125.3385008544858

length=40m

I tried to use this library

from geographiclib.geodesic import Geodesic

Geodesic.WGS84.Direct(start_point[1], start_point[0], bearing, length)

and my code

bearing=math.radians(bearing)
lat1=math.radians(start_point[1])
lon1=math.radians(start_point[0])
lat2=math.asin(math.sin(lat1)*math.cos(length/6378.1)+math.cos(lat1)*math.sin(length/6378.1)*math.cos(bearing))
lon2=lon1+math.atan2(math.sin(bearing)*math.sin(length/6378.1)*math.cos(lat1),math.cos(length/6378.1)-math.sin(lat1)*math.sin(lat2))
lat2=math.degrees(lat2)
lon2=math.degrees(lon2)

results are the same and both are incorrect
below is an example


texts:
How can I disable zoom after the locate button (widget locate) has successfully found user location?

I am using esri JS (arcgis version 3.20).
I tried looking at docs but could not find answer whatsoever.
https://developers.arcgis.com/javascript/3/jssamples/widget_locate.html


texts: Automatically populating ArcGIS Online map using ModelBuilder?

texts:
I have created a little model in ModelBuilder to take some cuts of some data e.g. last 30 days, last 6 months etc. I have then published this output to a feature service using ModelBuilder. I have set up all the names so it will overwrite itself meaning the model can just be ran daily and the new feature will be published and overwrite itself online, however it will not let me overwrite the service definition even though "overwrite processing" is turned on in geo-processing options. Any thought on why?

From this point I have to manually open up my web map and remove the old layer and add the new one into the map. Is their a way to incorporate doing this using my model as the reason I'm creating this is for any skill level user can run this model in my absence.

I am using ArcGIS 10.3 for Desktop but I've only got the basic license with the productivity suite


texts:
I am completely new to GIS, and am trying to extract data from PRISM Climate data.

I used the Point Sampling Tool in QGIS, but faced the following error:

  The new layer seems to be created, but is invalid. It won't be loaded

I have two layers:

My raster layer is precipitation data from PRISM.
My delimited text layer (.csv) has the following structure with a sample row:

city        state   latitude     longitude      source
marathon    texas   30.2066498  -103.2464676    Bing Maps

What could I be doing wrong here?

I previously tried sampling data from the same dataset (different time period), but with a different delimited text layer, and no such error occurred. Strangely, in my attempt to figure out this error, I tried it again and this error appeared.


texts:
http://imm.io/2cpJ

What is the thinking of OGC to keep this point that "point have no boundry , it just has interior"... what do they mean?

this notation: F**F*****...still i am confused due to OGC model.

It should be TFTTFTTFT as far as I  know: the point have no interior, it has just a boundary, so it should be an F and T relation between point and line.


texts:
Somehow I believe to recall there does exist some ind of "draw parallel line" plugin for QGIS.

Could anyone please point me to - I can't find it :-/


texts: Is there a source for free CSV files with geo data and other desired fields?

texts:
I'm writing some software components that process geospatial data. The input data is a CSV with general columns (e.g. name) and geo columns (e.g. city, lat or zip). I would like to test these components with some thousands of points.

Is there an online source that generates free csv files with geo data, according to a list of desired fields?

For example, I would like to have 1,000 points in North America, each with lat, lon, city, zip and state. Or 5,000 points in Europe, each with lat and lon.

I will probably build it myself from free data sources if it doesn't exist, but using an existing source would save lots of time.


texts: Downloading features from website with interactive map as shapefiles?

texts:
I found an interactive map on a website: here

The map on the website looks like this:

Is there any way to do download the red lines as shapefiles?


texts: Problem with coordinates not being calculated in decimal degrees

texts: I have a qGIS question regarding the calculation of x and y coordinates and showing the results in decimal degrees.  Initially I was using the 3.2 version, but went back to the 2.18.13 version because the newer version was really off for the coordinates (it would return 40.7768, -65.7466 when hovering over a point that was closer to 38.78, -76.87, for example).  The 2.18.13 version of qGIS seems to be giving the right coordinates when hovering over a point, but when calculating them into a new field, the results returned aren't in decimal degrees (i.e. 1343104.31085 for x and 531546.19274 for y).  The layer is set to a CRS of NAD83 HARN (EPSG:4269 (OTF)). Is there a conversion factor that should be used?  I have tried saving the dataset as a shp file (it was already), adding fields and calculating using $x and $y, and I have also used the toolbox to add geometry, and receive the same results - not in decimal degrees.  Any advice?


texts:
Does anyone know of a program, extension or tool that links a video with a GPS trail?

I found this video and would like to be able to do something similar.
http://www.youtube.com/watch?v=G07DWVRJBVU


texts: Limit Zoom on WMTS Layer using Openlayers 2.12

texts:
I'm implementing a baselayer to be used by a Plone addon (collective.geo.openlayers). My implementation specifically creates a layer and returns it to the addon, which enables many different layers to be configured through a webinterface.

I got it working well, but there is one problem I'm facing that I cannot wrap my head around.

The map displayed is a small local map which only shows tiles when zoomed in to zoom level 11 and higher. For this reason I would like to limit the user from zooming out farther than that.

But when I try to add maxZoomLevel, minZoomLevel or numZoomLevels to the layer definition I get type errors in the getTileInfo function because Openlayers is trying to access the longitude of the tileOrigin property which for some reason is null.

I tried different suggestions about limiting the zoom level but to no avail.

Does anyone know how to limit the zoom levels in the following layer definition?

I cannot simply work around the issue by doing something to the map instance itself, as I only control the layer creation.

function() {

var zugmap_matrixids = [{
        identifier : 0,
        topLeftCorner : new OpenLayers.LonLat(919690.324327, 5987771.047746)
    }, {
        identifier : 1,
        topLeftCorner : new OpenLayers.LonLat(929474.2639475, 5987771.047746)
    }, {
        identifier : 2,
        topLeftCorner : new OpenLayers.LonLat(929474.2639456, 5987771.04773376)
    }, {
        identifier : 3,
        topLeftCorner : new OpenLayers.LonLat(931920.24885072, 5985325.06282864)
    }, {
        identifier : 4,
        topLeftCorner : new OpenLayers.LonLat(931920.24885072, 5984102.07037608)
    }, {
        identifier : 5,
        topLeftCorner : new OpenLayers.LonLat(932531.745078525, 5983490.57415958)
    }, {
        identifier : 6,
        topLeftCorner : new OpenLayers.LonLat(932837.49319014, 5983184.82603666)
    }, {
        identifier : 7,
        topLeftCorner : new OpenLayers.LonLat(932990.36724671, 5983031.95198009)
    }
];

var layer = new OpenLayers.Layer.WMTS({
    name : "LuftbildPlus+",
    url : [
        "http://webdienste.zugmap.ch/OrthofotoPlus_tiled/service.svc/get"
    ],
    layer : "LuftbildPlusGoogle",
    matrixSet : "LuftbildPlusGoogle",
    matrixIds : zugmap_matrixids,
    tileFullExtent : new OpenLayers.Bounds(919690.324327, 5948635.289264, 978393.96205, 5987771.047746),
    maxExtent : new OpenLayers.Bounds(919690.324327, 5948635.289264, 978393, 5987771),
    style : "default",
    format : "image/png",
    isBaseLayer : true,
    zoomOffset : -11,
    //numZoomLevels: 8
});

return layer
}

You'll notice that numZoomLevels is commented out. I feel like this should be the correct property to set as the matrices only contain 8 levels, but if I uncomment it, I get the type error I was talking about. If I don't use numZoomLevels or any other property that controls the zoom the map works fine.


texts: Iterate variables in Con statement in Raster Calculator- model builder

texts:
I have this formula which works fine in raster calculator on its own. Here is a simplified version of it:

Con((("Trial_analysis\winter_pm" == 2)) & ("Wind_directions\winter_pm_WD" <=b) & ("Wind_directions\winter_pm_WD" >= c) ,d)

I need however to iterate through the variables (b,c,d), I have tables for all these variables. Basically I need the model to go through the winter_pm raster, iterate through all the conditions, and generate rasters from the list of values 'd'. So for every  b and c that are true, generate a raster with value d (of course limited to the cells where those conditions are true).
I tried implementing many variations in model builder using iterate value fields then calculate values for the tables a,b,c,d. But when I use these as inline variables in the raster calculator condition statement only ONE raster is generated. Any help would be great!


texts: ESRI Search Widget - Alter Search Parameters (ArcGIS JS v4.9)

texts: The Search widget within the ArcGIS JavaScript API (v 4.9) works by returning all records that contain your search string.  In other words, if you search for "123" You'll get "123 S Main St", "5123 W Oak St", "Painting 123", and more.  I just want to limit the search to only those records starting with the search string, not everything containing it.  There seems to be no exposed way however to alter what appears to be the "LIKE %search string%" query baked into the widget.  Has anyone been able to accomplish this?


texts: Is it possible to replicate an SDE geodatabase including the archive tables?

texts:
I have a database in one location which I would like to replicate to another using one-way replication. In testing we haven't been able to replicate the archive tables, which is unfortunate because we need to keep these historical versions in both locations.

Is it possible to simply replicate the whole DEFAULT database (including archive tables).

If not, is there any way to work around this?


texts: Auto-selecting records in all spatial and non-spatial joined/related tables by selecting only "parent" spatial record in ArcGIS for Desktop?

texts:
(1) I have a layer of pipes with related spatial and non-spatial data linked via joins and relates. Is there a way of selecting the top/parent record and then having ArcGIS auto-select all related records down the tree? Seems to work for spatially joined tables but not non-spatial ones.

(2) At a later stage I would also like to insert a step before this, of having a user prompt where someone can enter an asset ID and step (1) above will be done, with only the linked records shown in each table, where the end-user can then click on the different table tabs at the bottom of the attribute table window and only see the linked records.


texts: Adding GPX vector layer not loading all link elements in file

texts:
I am loading a GPX file into QGIS 2.18 and all seems to work ok. I am finding that the attributes table is only showing me the first two link_href entries. I have some points and tracks in my GPX file that have more than 2 link entries and these seem to be getting dropped on import.

Is there a way to extend the capabilities of the GPX import to load up more of the link entries in my GPX files? Adding 3 more to allow up to 5 links should be all that I need.


texts: Delete empty and null fields using arcpy

texts:
I need to remove fields within file GDB and SDE databases that contain only empty or null values. I was hoping to do this for an entire database at a time rather than per feature class. The code below shows the best of what I have so far. It does not, however, accomplish what I had hoped. I am somewhat new to Python and am using ArcMap 10.3.1.

import arcpy
from arcpy import env
env.workspace = r"C:\New Folder\Default.gdb"

def DelFields(path):
    env.workspace = path
    # first need to point to FCs within DSs
    datasets = arcpy.ListDatasets(feature_type='All')
    datasets = [''] + datasets if datasets is not None else []

    for ds in datasets:
        for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
            # list comp syntax:
            # [{code}var{code} for var in [list] if var{logic test}].
            fieldList = [f.name for f in arcpy.ListFields(fc) if f.type ==
                         "Short" or "Long" or "Float" or "Double" or "Text" or
                         "Date" or "Blob" or "Raster" or "Guid"] # "All" did not work
            print fieldList # to see what I have so far

            if fieldList: # this is a bit of a long shot I think
                # The UpdateCursor line gives error:
                # TypeError: 'field_names' must be string or non empty sequence of strings
                with arcpy.da.UpdateCursor(fc, [fieldList]) as cursor:
                    for row in cursor:
                        for i in range (len(fieldList)):
                            if row[i].rstrip() == "": # trying to find empties
                                arcpy.DeleteField_management(path, fieldList)
                            elif row[i] == "<Null>": # trying to find Nulls
                                arcpy.DeleteField_management(path, fieldList)
                                cursor.updateRow(row)

DelFields(r"C:\New Folder\Default.gdb")


texts: WMS + Google Earth + tilted view = alignment failure?

texts:
I'm working on a custom web-based application that generates spatial data and that includes a web-based map viewer.  We have a custom servlet that handles WMS requests (not a complete WMS implementation, just enough for our needs) and uses MapServer/SWIG to generate image data.  The web front end uses OpenLayers.  So far, so good: everything works and gives us nice-looking interactive maps.

We also built a feature that creates a KML file with a <GroundOverlay> that refers to our WMS URL -- the user clicks a link and it opens in Google Earth and they can see the data.  Again, so far, so good.

The problem comes when the user starts to tilt the camera and pan around in Google Earth.  Google Earth starts mis-registering our data, at first by just a few meters at a time but eventually by tens of kilometers, if you tilt, zoom in and out, and pan around enough.

I did some searches and I found several sites that indicated that Google Earth's support for tilting WMS imagery was poor, but little in the way of specifics.  Has anyone else tried this and experienced similar issues?

More importantly, how should I approach this problem?  Given that our users are very interested in doing exactly this kind of exploration (plot data with topography, tilt the camera down, and fly around in our subject area) and that our users are already pretty committed to using Google Earth as the client, how should we make the data available so that they can view it?


texts:
I recently uninstalled the ArcGIS 10 from my computer, the software got uninstalled , but I am unable to uninstall the License manager. It is giving me an error  "you don't have sufficient access to uninstall ArcGIS license manager 10"..and I am logged in as Admin.

Due to this problem I am not able to install the ArcGIS 10.1.


texts:
The Vector menu no longer works.  All the other menus drop down and show the different options available. Here is the error message shown upon startup. I do not know if this is linked to the Vector menu or if two different problems.

2019-03-24T14:06:40     WARNING    Python error : Couldn't load plugin 'NITK_RS-GIS_17' due to an error when calling its classFactory() method See message log (Python Error) for more details.
2019-03-24T14:06:42     WARNING    Plugin processing : The plugin will be disabled because it crashed QGIS during last startup. Please report an issue and re-enable the plugin when the problem has been solved.


texts: How to automatically get mode (modus) of a raster in ArcGIS

texts:
do you have some experience how to automatically get mode value from a raster ( value that appears most often in a dataset)? I am looking for something like "get raster properties" in Spatial analysis in ArcGIS, unfortunately there is no possibility to directly calculate mode. Subsequently, I would like to use this calculation of mode in Model Builder. Based on "mode value" I would like to reclassify raster image into binary file
(0= values < mode, 1 = values > mode).


texts: Label not inside border after map rotation

texts: Does anyone know why my label doesn't show normally within the label border when I apply a rotation to my map? It shows fine with zero rotation.


texts:
I have created two caches for a project and there are coexisting within the same folder.

Right behind I open ArcMap, add my File System Raster Layers and Share as Service.
My *.mxd document has exactly the name as the hosting file from my caches folder.
ArcGIS Server detect automatically the caches and I can publish the service.
But after that operation done within manager the cache seems to be not consume.

Any idea how to make it happen?


texts: Extracting unique values from raster in area using ArcGIS Desktop?

texts:
I have Landcover raster and thousand points. I want to create buffer  around each point a extract all unqiue values in buffer from Landcover raster.

Do you know how can it be done using ArcGIS Desktop?

An example of what I want to create is in this picture.


texts:
QGIS stopped working with this error:

/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libgdal.so.20: undefined symbol: GEOSDelaunayTriangulation_r

I removed/purged everything for a fresh install but I continue to get the same error, I have tried multiple different sources. I can successfully install QGIS 2.8 but 2.14 and 2.18 do not work.

This is the output of attempting to install it with aptitude.

 Setting up qgis-providers (2.14.20+dfsg-1~xenial0) ...
/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libgdal.so.20: undefined symbol: GEOSDelaunayTriangulation_r
dpkg: error processing package qgis-providers (--configure):
 subprocess installed post-installation script returned error exit status 127
dpkg: dependency problems prevent configuration of qgis:
 qgis depends on qgis-providers (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis-providers is not configured yet.

dpkg: error processing package qgis (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-provider-grass:
 qgis-provider-grass depends on qgis (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-provider-grass (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-plugin-grass:
 qgis-plugin-grass depends on qgis (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis is not configured yet.
 qgis-plugin-grass depends on qgis-provider-grass (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis-provider-grass is not configured yet.

dpkg: error processing package qgis-plugin-grass (--configure):
 dependency problems - leaving unconfigured
Errors were encountered while processing:
 qgis-providers
 qgis
 qgis-provider-grass
 qgis-plugin-grass
E: Sub-process /usr/bin/dpkg returned an error code (1)
Failed to perform requested operation on package.  Trying to recover:
Setting up qgis-providers (2.14.20+dfsg-1~xenial0) ...
/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libgdal.so.20: undefined symbol: GEOSDelaunayTriangulation_r
dpkg: error processing package qgis-providers (--configure):
 subprocess installed post-installation script returned error exit status 127
dpkg: dependency problems prevent configuration of qgis:
 qgis depends on qgis-providers (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis-providers is not configured yet.

dpkg: error processing package qgis (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-plugin-grass:
 qgis-plugin-grass depends on qgis (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-plugin-grass (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-provider-grass:
 qgis-provider-grass depends on qgis (= 2.14.20+dfsg-1~xenial0); however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-provider-grass (--configure):
 dependency problems - leaving unconfigured
Errors were encountered while processing:
 qgis-providers
 qgis
 qgis-plugin-grass
 qgis-provider-grass

I must've updated some part of my system that is causing gdal? to fail.
I am on Linux Mint 18.1


texts: Get shortest/fastest path between two points from SpatialLines(DataFrame)

texts:
I have a large rivernetwork as a SpatialLines(DataFrame) and want to subset it so I get a SpatialLine with only the shortest/fastest path between two SpatialPoints(DataFrame) which propably need to be snapped onto the line first.

In the example Code, I want to get the closest line between the red dots, while moving on the lines. Also I want to move between the left red point and the right blue point, which needs to be snapped to the closest SpatialLine

library(sp)
x <- c(1,5,6,8)
y1 <- c(1,3,4,7)
y2 <- c(5,5,5,2)

L <- SpatialLines(list(Lines(Line(cbind(x,y1)), ID="a"),Lines(Line(cbind(x,y2)), ID="b")))
P <- SpatialPoints(data.frame(x=c(1,8),y=c(1,2)))
P_snap <- SpatialPoints(data.frame(x=c(8),y=c(1)))

plot(L)
points(P,col="red")
points(P_snap,col="blue")


texts:
I'm totally new to GIS.  However, I've proudly managed to create a map, which is correctly geo-referenced with a Great Britain Ordnance Survey (OS) background (streetview level) and I've added some point data.

I would like to know how to create or import a layer with the OS grid lines.  I'm sure this has been done a million times, but I can't seem to find the answer on the web - I suspect I'm not using the correct search terms.

Any help appreciated.


texts: How to apply categorized symbols to polygon layer with PyQGIS?

texts:
I'm fetching a polygon layer from postgres using the following:

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "mydb", "postgres", "xxx")
uri.setDataSource("public", "my_table", "geom")
lyr = QgsVectorLayer(uri.uri(), "Blocks", "postgres")
QgsMapLayerRegistry.instance().addMapLayer(lyr)

This layer has an attribute field named blocks that has 6 values:
    'A-6', 'B-4', 'G-5', 'D-2', 'P-7', 'C-3'.

I want to apply the color green to polygons where the attribute 'block' value is 'B-4' and red color to the other polygons in the layer that have a different value in the 'block' attribute.

I've seen code to apply graduated symbols using PyQGIS but I want categorized symbols.

How can I do this?


texts: Grouping points into user-provided number of clusters using ArcGIS for Desktop?

texts: I am looking for a process in ArcGIS that takes a point feature class and groups the points into a user-provided number of clusters.


texts: Auto-updating count of points in polygons?

texts: I have a Point shapefile and a Polygon shapefile. Is it possible to create a field into the attribute table of the Polygon shapefile with the count of points inside each polygon, such that if I add a new point inside a polygon in a later moment the change is automatically reflected in the Count Field?


texts: How to repeat drawing function in Leaflet

texts:
I am trying to create a Leaflet.js web map that gives users the opportunity to draw lines.
this line will be sent to mapquest elevation API to be retun to the map as a GeoJSON with elevation data calculated in order to showing some info about elevation in a popup for that drawn line.

The function goes well, but when i try to draw another line it doesn't work well and returns the first drawn line.

Here is my code:

    map.on(L.Draw.Event.CREATED, function(event) {
    var layer = event.layer;
    // var content = getPopupContent(layer);
    // if (content !== null) {
    //     layer.bindPopup(content);
    // }
    drawnItems.addLayer(layer);
    drawData = drawnItems.toGeoJSON();

        function showElevation(layer) {
        var coords = layer.feature.geometry.coordinates;
        var minPoint = minHPoint(coords, (layer.feature.geometry.type == 'LineString'));
        var maxPoint = maxHPoint(coords, (layer.feature.geometry.type == 'LineString'));
        minH = minPoint[2];
        maxH = maxPoint[2]

      if (elevationDisplayed && (currFeature === layer.feature)) {
        map.removeControl(el);
        elevationDisplayed = false;
        }
      else {
        el.clear();
        if (!elevationDisplayed) el.addTo(map);
        el.addData(layer.feature);
        currFeature = layer.feature;
        elevationDisplayed = true;
      }
    };
        var profile = drawData.features[0].geometry.coordinates;
        fetch('http://open.mapquestapi.com/elevation/v1/profile?key=tHXSNAGXRx6LAoBNdgjjLycOhGqJalg7&shapeFormat=raw&latLngCollection='+profile)
          .then(r => r.json())
          .then(data => {
            var latlngs = [];
                for (var i=0; i<data.shapePoints.length; i=i+2) {
                    latlngs.push([data.shapePoints[i], data.shapePoints[i+1], data.elevationProfile[i/2].height]);
                }

            geojson = {
                "type": "Feature",
                "geometry": {
                    "type": "LineString",
                    "coordinates": latlngs
                }
            };

            var popupGeojson = L.popup();

            L.geoJSON(geojson, {
                style: function(feature, layer) {
                    return {'color': 'red'}
                },
                onEachFeature: function(feature, layer){
                    showElevation(layer);
                    $('#hideElevation').show();
                    $('#hideElevation').click(function() {
                        $('.background').toggle();
                    });
                    var firstPoint = feature.geometry.coordinates[0][2];
                        gainElevation = maxH - firstPoint;
                        var length = turf.length(geojson, {units: 'kilometers'});
                        popupGeojson.setContent("Distance: "+''+Math.round(length)+ " km"+''+'
'+"Low point: "+''+minH+' meters'+''+'
'+"High point: "+''+maxH+' meters'+''+'
'+"Gain: "+''+Math.round(gainElevation)+" meters"+'');
                        //layer.bindPopup(popupGeojson);
                        map.on("zoomend", function() {
                            layer.bindPopup(popupGeojson).openPopup();
                        });
                        map.fitBounds(layer.getBounds());
                }
            }).addTo(drawnLine);
            drawnLine.addTo(map);
          })
});