texts: Iterating over memory layer to run processing algorithm using Python in QGIS?

texts:
I am trying to write a python script that iterates over the features of a layer, selects a features and feeds into processing.runalg(). I am pretty sure the selection is going on. However, the processing.runalg("qgis:xxx") bit ocurrs oblivious to the selection! I know there are plenty of similar questions and answers, but I have tried their suggestions and haven't been successful. I though maybe if someone sees the code maybe could help me resolve this.

Here is the code (indents in the for loops might be messed up from formatting):

inter_shpmemory = iface.addVectorLayer(main_directory + "\\" +
intersect_shpname,"Zone layer",labtype)
specieslayer_mem = iface.addVectorLayer(main_directory + "\\" +
bd_shpname,"GBIF layer",labtype)

###Create a folder to safe all stuff
newfolder = main_directory + "\\" + "Zonal_GBIF"

if not os.path.exists(newfolder):
    os.makedirs(newfolder)

it = specieslayer_mem.getFeatures( QgsFeatureRequest() )
inputLayerIds = [i.id() for i in it]
#selects all of the input feature and one feature from the clip layer
specieslayer_mem.setSelectedFeatures( inputLayerIds )###is this even
#needed?

features = inter_shpmemory.getFeatures()

for feature in features:
            indbd_name = feature['simp_name']
            inter_shpmemory.setSelectedFeatures([feature.id()])
            print "Feature ID %d: " % feature.id()
            print indbd_name
            processing.runalg('qgis:clip', specieslayer_mem, inter_shpmemory, newfolder + '\\' + indbd_name +'.shp')
            inter_shpmemory.setSelectedFeatures([])

specieslayer_mem.setSelectedFeatures([])

I can't get to do it what I want.


texts:
I was able to run my application on using Legend earlier but now I want to update my Legend when scrolling on the map. So I tried to used this code GetLegendGraphic and getLegendUrl function filter by extent . And I got the gist of how he is able to do that but when I am trying the same I keeping getting the error.
var graphicUrl = wmsSource.getLegendUrl(resolution);

This is mainly due to the getLegendUrl is not getting declare. I have tried to initialize the same using different ways but it seems getLegendUrl is not getting declare correctly.
This is my JS
var wmsSource = new ol.source.TileWMS({
    url: 'http://localhost:8089/geoserver/India/wms',
    params: {'LAYERS': 'India:India', 'TILED': true,'TRANSPARENT':true},
    ratio: 1,
    serverType: 'geoserver',
    transition: 0,
})

var getExtent = function () {
    return map.getView().calculateExtent();
}
var getProjectionCode = function () {
    return map.getView().getProjection().getCode();
}

var updateLegend = function(resolution) {
    console.log(wmsSource);
    console.log(wmsSource.getLegendUrl(resolution));
    var graphicUrl = wmsSource.getLegendUrl(resolution);

    const crs = `&CRS=${getProjectionCode()}`;
    const bbox = `&BBOX=${getExtent().join(',')}`;
    const height = getHeight(getExtent());
    const width = getWidth(getExtent());
    const heightAndWidth = `&srcwidth=${height}&srcheight=${width}` ;
    const legenOptions = "&legend_options=countMatched:true;fontAntiAliasing:true;hideEmptyRules:true;forceLabels:on" ;

    var img = document.getElementById('legend');
    img.src = `${graphicUrl}${crs}${bbox}${heightAndWidth}${legenOptions}`;
};

var map = new ol.Map({
      target : 'mymap',
      view: myview,
      layers: [
        new ol.layer.Tile({
            source : new ol.source.OSM({
                attributions: 'Surveyor Allpication'
            })
        }),
        new ol.layer.Tile({
            source : wmsSource
        })
      ]
})

var resolution = map.getView().getResolution();
updateLegend(resolution);

This is my HTML

        Legend:

Any idea how can I make this work?
I keep getting the error of: wmsSource.getLegendUrl is not a function at updateLegend


texts:
I have compiled QGIS from branch QGIS-release-3.2 on bionic beaver (ubuntu 18.04), the ccmake step detects python 3.6 path (/usr/bin/python3.6) and all the compilation goes OK

Starting QGIS 3.2.2 the following error pops up

Reporting that it can't import traceback, more interesting is that there is no indication of python interpreter and/or paths being used

Inside the python terminal I can do  import qgis  without problems.

Question: How can I fix this could it be that somewhere QGIS doesnt have python path properly set ???

Note: The OS (bionic beaver) was upgaded from xenial and sometimes strange things happen


texts: I am working with a DEM of the lunar south pole in QGIS. I have currently set the projection to "no projection", or Cartesian co-ordinates (origin at the pole). The horizontal locations of features are consistent with other sources, but all of the elevation data appears to be exactly twice what it should be. Is there a way to re-scale the elevation data only, or are my settings in QGIS wrong?


texts: Running tool over group of layers using ArcGIS Pro

texts:
I want to run basic tools (e.g., Dissolve, Intersect) over a group of feature classes from the content pane.
Is there any way to have a tool run over a set of layers grouped within the contents pane?
I'm currently doing it layer by layer within the tool dialogue box. Even with Python, it seems I need to type out the names of each layer which can take longer than simply clicking from the drop-down menu in the toolbox.


texts:
I have an featureCollection Object and I want to fill it with values through a loop that takes values from another object. Here is what I did:

  for (var key in data){
                 var geoJsonFeatureCollection = {
                     type: 'FeatureCollection',
                     features:
                     {
                         type: 'Feature',
                         geometry: {
                             type: 'Point',
                             coordinates: [data[key].YO, data[key].XO]
                         },
                         properties: {
                             "origin_city": data[key].Origin,
                             "origin_lat": data[key].YO,
                             "origin_lon": data[key].XO,
                             "destination_city": data[key].Destination,
                             "destination_lat": data[key].Y,
                             "destination_lon": data[key].X

                         }
                     }};
                 }

But I need to place the loop inside the geoJsonFeatureCollection, how to do this ?


texts: QGIS-Python Need to check if layer is in current extent and make it active

texts:
I'm putting together a plugin for QGIS. The user inputs a county and a coordinate pair. The plugin then adds in the county data and zooms to the given coordinates.

The problem is that there are multiple counties with the same name, when zoomed into the coord pair, the county data may be visible but not active because there are multiple layers with the same name added.

I need a way to activate the layer that is in view/visible/can-be-seen-at-the current-position in canvas. I've tried:

import qgis
import qgis.utils
canvas = qgis.utils.iface.mapCanvas()

layers = qgis.utils.iface.legendInterface().layers()
i = qgis.utils.iface
legend = i.legendInterface()

#iter through the legend layers
for layer in layers:
#if the name of the layer is the same as the what the user input
if layer.name().lower() == county.lower():
    #and you can see it in the canvas/in the extent/if youre looking at it
    if canvas.extent() == layer.extent():
        #make the layer active so i can select, or identify things
        iface.setActiveLayer(layer)

legend.isLayerVisible(layer) is not what im looking for, it toggles layers, the data is already visible, I need the layer in my current position active.


texts: OpenLayers SelectFeature control ignores graphicZIndex

texts:
I am having trouble with SelectFeature control with hover = true and highlightOnly = true.

There are two adjacent polygons on the same layer, one on top of the other. There's also a SelectFeature control that highlights a polygon on hover. Layer has styleMap with 'default' and 'select' intents. Select intent has a higher graphicZIndex. But it's ignored - whenever I hover upper polygon its lower border (the one that's common with the lower polygon) is rendered under the lower polygon.

If you go to http://jsfiddle.net/vLpt6/12/ you'll see two adjacent polygons - hover the upper polygon and its lower border will be rendered under the lower polygon as if graphicZIndex had no effect.

Here's the HTML code:

And javascript code:

var map = new OpenLayers.Map({
    div: 'map',
    allOverlays: true,
    controls: [],
    fractionalZoom: true
});

var layer = new OpenLayers.Layer.Vector('Layer');
map.addLayer(layer);

layer.rendererOptions = { zIndexing: true };
layer.styleMap = new OpenLayers.StyleMap({
    'default': {
        fillColor: '#a9a8a7',
        strokeColor: '#ffffff',
        strokeWidth: 1,
        graphicZIndex: 10
    },
    'select': {
        fillColor: '#a9a8a7',
        strokeColor: '#5d7cc4',
        strokeWidth: 5,
        cursor: 'pointer',
        graphicZIndex: 100
    }
});

var format = new OpenLayers.Format.WKT();
var features = [
    format.read('POLYGON((0 0, 100 0, 100 50, 0 50, 0 0))'),
    format.read('POLYGON((0 0, 0 -50, 100 -50, 100 0, 0 0))')
];
layer.addFeatures(features);

map.zoomToExtent(layer.getDataExtent());

var hoverControl = new OpenLayers.Control.SelectFeature(
    [layer],
    {
        hover: true,
        highlightOnly: true
    }
);
map.addControl(hoverControl);
hoverControl.activate();

What should I correct so that graphicsZIndex would have effect?map


texts: Combining hyperspectral bands to form a multispectral image?

texts:
I understand that the difference between the two are the number and width of the bands. Given a hyperspectral image, is it possible to 'convert' it to multispectral ?

In this paper for example, the authors make use of RGB + Near Infrared (multispectral) for scene recognition. If i have a hyperspectral image and I keep only 4 bands (3 bands for rgb and 1 for NIR) while discarding/ignoring the rest, would I be able to achieve a somewhat similar result or would too much information be lost since the bands are much narrower than multispectral bands ?

Or should I perform some algorithm to combine the bands from the NIR and RGB respectively such that I have a proper multispectral image ?


texts:
In QGis 2.18.3

The Vector -> Geometry Tools -> Extract Nodes tool will generate a point layer shapefile with points at each vertex of each line.

But each point in layer does not have the z coordinate, only x and y.

How can i extract nodes with (X Y Z) from lines?


texts: I have two Landsat 8 satellite images from two different dates in may. When I import them in ArcGIS they slightly overlap but the pixel colours don't match making it hard to run an NDVI. I have been trying to run colour matching on the two datasets by creating a mosaic or raster catalogue but it doesn't work.


texts:
I have a shapefile that I was working on in ArcGIS Desktop 9.3 get corrupted somehow and it gives me the error "Number of shapes does not match the number of table records.". What I would want is to recover this shapefile as it has a lot of edits that other revisions of the file do not have.

Earlier I recovered from this error using ogr2gui but the tool now just crashes. I have even tried the command line version ogr2ogr but still no luck. This script, Shapefile Repair Tool from the ArcScripts website didn't help either.

Any ideas? (except that I stop using shapefiles)

Here is the output from ogr2ogr (based on a a suggestion in the answers):


texts: How to bring LSOA data up to MSOA level London using ArcGIS for Desktop?

texts: I need to preform regression analysis but my issue is that I have data at 2 different Output areas, MSOA and LSOA. My question is, how to I bring my LSOA data up to MSOA level?


texts:
I need a way of comparing two tracks (they will most likely be GPX files), i.e., 2 arrays of lat-longs.  I need to write some code that will compare the two arrays and give me a measure of how close they are to one another.  I am testing various methods of determining location against a known route and would like to compare them.

Would the sum of the cross-track errors be sufficient?

I have seen some articles on Hausdorff distance (and dimension) and in particular Wikipedia describes "A measure for the dissimilarity of two shapes is given by Hausdorff distance up to isometry", which looks promising but I confess I don't understand what this means let alone how to implement the algorithm in code.

In short I would like to be able to state that method 1 is twice as good as method 2 at determining location.


texts: Why is ArcMap reverting 4-band raster settings?

texts:
ArcMap keeps redrawing my 4 band false colour image in standard 3 colour bands.  As you can see below I have my 4 band data set to auto display as (4,1,2 - infrared, red, green).  When I bring in an image layer it displays as it should and the layer properties clearly show my settings are working.

If however, I pan around and zoom in the image reverts to a 3 band image - the image difference is shown in the attached image.  When I go back into image properties it shows that ArcMap has automatically changed my bands to (1,2,3 - R,G,B) - even though my customization menu still shows my desired settings as (4,1,2 - iR, R, G).

Any idea why ArcMap is doing this?

I


texts: All layers in ArcMap lined up, but when exported to CAD, they shifted to different locations

texts:
I have all the layers in ArcMap lined up correctly, but they shifted in CAD after I exported them. Why is it?

The data frame property shows current coordinate system:

Current coordinate system:

  NAD_1983_HARN_StatePlan_Oregon_North_FIPS_3601_Feet_Intl

  Projection: Lamber_Conformal_Conic. Geographic Coordinate System:

  GCS_North_American_1983_HARN Datum : D_North_American_1983_HARN

The layer of the drainage system feature class shows:

Coordinate System:

  NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601_Feet_Intl

  Projection : Lambert_Conformal_Conic

  Geographic Coordinate System : GCS_North_American_1983_HARN Datum : D_North_American_1983_HARN

The contour layer feature class has a Projected Coordinate System:

  NAD_1983_2011_Oregon_Statewide_Lambert_Feet_Intl

  Projection: Lambert_Conformal_Conic

  Geographic Coordinate System : GCS_NAD_1983_2011 Datum :
  D_NAD_1983_2011


texts: How to load GeoTIFF info by gdal library like Matlab geotiffinfo function

texts:
geotiffinfo function in Matlab show different geotiff info from GeoTIFF images.

I use gdal library for read images info.

How to load GeoTIFF info by gdal2 like Matlab geotiffinfo function?


texts:
I am using google maps and using markers on the map. When user clicks on any marker showing related information in the Infowindow and every thing is working fine.

But every time the infowindow opens on top of the marker. Our requirement is infowindow should open/come left,right,top or bottom with respective to the available space.

Thanks

 google.maps.event.addListener(marker, 'click', function() {
         if(infowindow){
             infowindow.close();
         }
         infowindow  = new google.maps.InfoWindow();
         //var centerPos = new google.maps.LatLng(latlng.lat(), latlng.lng());
         //map.setCenter(centerPos,65);
         infowindow.setContent(contentString);
         infowindow.open(map,marker);

    });


texts: How can a get OSM data in QGIS in a way that shows the direction of a road as it would in JOSM?

texts: I am trying to download OSM data into QGIS in a way that is as close as possible to what it looks like in JOSM. I am doing a visual comparison between OSM data and satellite images. One of the things I am looking at is comparing if one-ways are mapped correctly. How can I get the OSM data in QGIS in a way that shows the direction of a road as it would in JOSM?


texts:
I am using EA DTM squares as input layers into the profile tool.

I have created a polyline and am using this as selection layer. A profile shape is correctly created.  The X axis is in meters and shows the scale correctly being the length of distance along the polyline.

However, the Y axis comes out from zero to 50000. In real life the elevation should be between 0 and 54 meters.

Any idea what the Y scale is being used and how this can be altered to show meters like the X axis? Both the DTM rasters and Polyline are using OSGB 1936 / British National Grid.


texts: Coupled Hidden Markov Model Matlab Error

texts:
I want to train CHMM using Matlab and I am new to use this software. I am using BNET CHMM which is this:

function bnet = mk_chmm(N, Q, Y, discrete_obs, coupled, CPD)
% MK_CHMM Make a coupled Hidden Markov Model
%
% There are N hidden nodes, each connected to itself and its two nearest neighbors in the next
% slice (apart from the edges, where there is 1 nearest neighbor).
%
% Example: If N = 3, the hidden backbone is as follows, where all arrows point to the right
%
% X1--X2
%   \/
%   /\
% X2--X2
%   \/
%   /\
% X3--X3
%
% Each hidden node has a "private" observed child (not shown).
%
% BNET = MK_CHMM(N, Q, Y)
% Each hidden node is discrete and has Q values.
% Each observed node is a Gaussian vector of length Y.
%
% BNET = MK_CHMM(N, Q, Y, DISCRETE_OBS)
% If discrete_obs = 1, the observations are discrete (values in {1, .., Y}).
%
% BNET = MK_CHMM(N, Q, Y, DISCRETE_OBS, COUPLED)
% If coupled = 0, the chains are not coupled, i.e., we make N parallel HMMs.
%
% BNET = MK_CHMM(N, Q, Y, DISCRETE_OBS, COUPLED, CPDs)
% means use the specified CPD structures instead of creating random params.
%  CPD{i}.CPT, i=1:N specifies the prior
%  CPD{i}.CPT, i=2N+1:3N specifies the transition model
%  CPD{i}.mean, CPD{i}.cov, i=N+1:2N specifies the observation model if Gaussian
%  CPD{i}.CPT, i=N+1:2N if discrete

if nargin < 2, Q = 2; end
if nargin < 3, Y = 1; end
if nargin < 4, discrete_obs = 0; end
if nargin < 5, coupled = 1; end
if nargin < 6, rnd = 1; else rnd = 0; end

ss = N*2;
hnodes = 1:N;
onodes = (1:N)+N;

intra = zeros(ss);
for i=1:N
  intra(hnodes(i), onodes(i))=1;
end

inter = zeros(ss);
if coupled
  for i=1:N
    inter(i, max(i-1,1):min(i+1,N))=1;
  end
else
  inter(1:N, 1:N) = eye(N);
end

ns = [Q*ones(1,N) Y*ones(1,N)];

eclass1 = [hnodes onodes];
eclass2 = [hnodes+ss onodes];
if discrete_obs
  dnodes = 1:ss;
else
  dnodes = hnodes;
end
bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'eclass1', eclass1, 'eclass2', eclass2, ...
          'observed', onodes);

if rnd
  for i=hnodes(:)'
    bnet.CPD{i} = tabular_CPD(bnet, i);
  end
  for i=onodes(:)'
    if discrete_obs
      bnet.CPD{i} = tabular_CPD(bnet, i);
    else
      bnet.CPD{i} = gaussian_CPD(bnet, i);
    end
  end
  for i=hnodes(:)'+ss
    bnet.CPD{i} = tabular_CPD(bnet, i);
  end
else
  for i=hnodes(:)'
    bnet.CPD{i} = tabular_CPD(bnet, i, CPD{i}.CPT);
  end
  for i=onodes(:)'
    if discrete_obs
      bnet.CPD{i} = tabular_CPD(bnet, i, CPD{i}.CPT);
    else
      bnet.CPD{i} = gaussian_CPD(bnet, i, CPD{i}.mean, CPD{i}.cov);
    end
  end
  for i=hnodes(:)'+ss
    bnet.CPD{i} = tabular_CPD(bnet, i, CPD{i}.CPT);
  end
end

I put the parameters in mk_chmm(4, 4, 4, 1, 1, 0) in the editor window, but it gives me an error:

Undefined function 'mk_dbn' for input
arguments of type 'double'.

Error in mk_chmm (line 70)
bnet = mk_dbn(intra, inter, ns,
'discrete', dnodes, 'eclass1', eclass1,
'eclass2', eclass2, ..

Any ideas how can I solve it please?? Sorry I am a beginner.


texts: Projection not changing in ArcMap after Project tool used?

texts:
I tried changing the shapefiles by projecting them by only using the Project tool and I made sure that the data frame properties were set to the correct projection.  I see the properties of the shapefile and the data frame properties say that they are all set to Albers Equal Area Conic but when they are displayed in a map it doesn't look like it.  Here is the image that I see.

Here is what one of my three source tabs looks like.  They all have the same projected coordinate system and even the same geographic coordinate system.


texts: PyQGIS categorized symbol renderer sort by label

texts:
I wan't to sort a QGIS categorized symbol renderer by label name with PyQGIS. First i create categories and after i apply them to my layer but how to make an ascending alphabetic order on categorie's name ?
Documentation explain we can sort categories by their label with .sortByLabel().
I'm using QGIS3.10.2.
# Categories
categories = []
for myAttrName, (color, label, contour) in myStyleDict.items():
    symbol = QgsFillSymbol.defaultSymbol(myLayer.geometryType())
    symbol.setColor(QColor(color))
    symbol.symbolLayer(0).setStrokeWidth(0.26)
    symbol.symbolLayer(0).setStrokeColor(QColor(contour))
    category = QgsRendererCategory(myAttrName, symbol, label)
    categories.append(category)

# Renderer
expression = 'myAttrField'
renderer = QgsCategorizedSymbolRenderer(expression, categories)
#renderer = QgsCategorizedSymbolRenderer(expression, categories).sortByLabel() ?
myLayer.setRenderer(renderer)


texts: Looping through project bookmarks using PyQGIS

texts: How can I access a project's bookmarks using PyQGIS?


texts: Automatically Remove Duplicate Vertices in a Polyline

texts: I am trying to automatically identify identical vertices within a polyline and remove them without converting the vertices to points first.  I have converted the vertices to XY points and identified a significant number of identical vertices, greater than 4000, that will need to be removed.  I believe I can use itertools combined with a search cursor and an update cursor using the SHAPE@ variable to identify identical vertices and delete one.  I have tried various bits of code and methods, but with little success.  Using itertools "overlap" did not identify the identical points in XY so I believe I will need to use the itertools "equal", but am really a bit lost at this point. I am using ArcGIS 10.2.2 and do not have access to any extensions. I would like to incorporate this into a larger Python script I am writing.


texts: How to fix missing ippj-5.3.dll error in QGIS?

texts:
I try to install OSGeo4W but when I tried to run it appeared an error message about ippj-5.3.dll is missing...
does anyone know what happen?


texts: How to identify unknown coordinate system?

texts:
Now I'm trying to convert unknown coordinate system to Google Maps (WGS84?). How can I identify what kind of projection am I dealing with? I'm inspecting www.travelinesoutheast.org.uk. Route response returns json with:

name:"Epsom, Epsom Railway Station (Stop N)"
coords:"-29903.75584,5319390.31291"

name:"Epsom, Epsom Race Course (W-bound)"
coords:"-28583.17541,5323006.72201"

link to route search result


texts: ogr2ogr: don't import properties on each feature?

texts:
I have a 2MB shapefile that when converted to GeoJSON becomes 10MB:

$ ogr2ogr -f "GeoJSON" -t_srs EPSG:900913
  landfill.geojson Authorised_Landfill_Sites_010k.shp

I can reduce the size of the GeoJSON a bit by adding -lco COORDINATE_PRECISION=4 to the command.

However, I've noticed that there are an awful lot of properties in the GeoJSON that I don't need, so I'd like to remove them completely:

 { "type": "Feature", "properties": { "LIC_ADMIN": null, "LIC_NMBR": null, "LIC_IPPCR": "HP3532LV", "LIC_WML": 0.000000, "CUST_NMBR": 0.000000, "STATUS": "Ef     fective", "LIC_LTYPE": "5.2 A(1) a)", "LIC_NAME": "Yorwaste Limited", "LIC_SITE": "HAREWOOD WHIN LANDFILL", "SITE_NAME": "HAREWOOD WHIN LANDFILL", "SITE_BUI     LD": null, "SITE_STRT": "Harewood Whin Landfill", "SITE_AREA": "Tinker Lane", "SITE_TOWN": null, "SITE_CNTY": "Rufforth", "SITE_PCODE": "YO23 3RR", "TYPE_DE     SC": "WASTE LANDFILLING; >10 T\/D WITH CAPACITY >25,000T EXCLUDING INERT WASTE", "NGR": null, "T_REF": "NE4731", "CTROID_X": 453600.000000, "CTROID_Y": 4513     00.000000, "REGION": "North East", "AREA": "Yorkshire", "DATE_ISSUE": "2007\/06\/21", "LIC_EPR": "EA\\EPR\\BK0507IB" }, "geometry": { "type": "Polygon", "co     ordinates": [ [ [ -131830.27, 7161699.84 ]

Can I use ogr2ogr to remove these properties at the time of conversion, and just keep the coordinates, or should I post-process the GeoJSON to do this?


texts: How to add column length in a pgRouting PostGIS database?

texts: I have 50 3d linestrings in postgis database, and I need to add column length so I could use pgRouting for calculating shortest path with Dijkstra. How can I add it?


texts:
I am wanting to utilize a high accuracy GNSS position in a GIS data collection application optimized for Chrome running on my Android device. The problem is that the device doesn't have Google Play or Google Play services, which seems to restrict the ability for Chrome to receive the GPS data from the location manager.

Has anyone experienced this in the past? If so, is there any work around available so Chrome can get the GPS position from the location manager?


texts: arcpy.mapping.ExportToTIFF fails to render layer?

texts:
Has anyone else had this problem?

I have an mxd file with two layers in it, one of which is fairly 'heavy' (i.e. takes a minute or so to draw). All layers are visible by default.

I have the following code, which opens the mxd file (from outside arcmap) and exports it to a tiff.

import arcpy
template = r"C:\mxd_templates\template1.mxd"
mxd = arcpy.mapping.MapDocument(template)
arcpy.mapping.ExportToTIFF(mxd, "C:/temp/test.tif")

The tif file is created, but the second 'heavy' layer is never rendered.
This happens no matter what image format I export to.

I've not found anywhere to be able to report bugs online, and I'm having a hard time to get ESRI to recognise this as a bug or to point out what I am doing wrong. Although it is hard to see what could possibly be wrong in 4 lines of code?

Same behaviour on 10.1, 10.2 and 10.2.2, 64bit. The layer that does not render is a raster in arcSDE.
The layer that does render is a raster in a file GDB

In response to comments below, Yes the map does export correctly when rendered manually via arcmap, and yes the layers are correctly ordered and have been thoroughly checked.


texts:
I added a function to display a popup on hover to the feature layer but it doesn't work very well.

Here my code:

map.on('pointermove', function(evt) {
  if (evt.dragging) {
    return;
  }

  var info = $('#info');
  info.html("");
  var pixel = map.getEventPixel(evt.originalEvent);
  var feature = map.forEachLayerAtPixel(pixel, function(feature) { return true; }, null, function(layer) { return layer === pmfeatlayer; });
  var viewResolution = map.getView().getResolution();
  var coordinate = evt.coordinate;
  url = pmfeatlayer.getSource().getGetFeatureInfoUrl(coordinate, viewResolution, projection, {'INFO_FORMAT': 'text/html'});

  if (feature && url) {

    info.css({
      left: pixel[0] + 'px',
      top: pixel[1] + 'px'
    });
//    setTimeout(function() {
    var xhttp = new XMLHttpRequest();
    xhttp.open("GET", url, true);
    xhttp.send();
    xhttp.onreadystatechange = function (aEvt) {
      if (xhttp.readyState == 4 && xhttp.status == 200) {

          var parser = new window.DOMParser();
          var res = parser.parseFromString(xhttp.responseText, "text/xml");
          var myitems = res.getElementsByTagName('item');

          if(res.getElementsByTagName('item').length != 0){
            info.html("
" + hover(xhttp.responseText) + "
").delay( 1000 ).fadeIn( 400 );
          } else {
            info.html("");
          }

       } // end if 4 && 200
     }; // end onreadystatechange
 //}, 1000);
   } // end if (feature && url)
 }); // end pointermove

As you can see I'm using map.on('pointermove'... But it's not the exact event, What can I use to detect that the cursor is above the feature layer and it's not moving, for example, .5 seconds?

I tried to add a delay function but still, it doesn't work properly, it made a lot of requests and the popup display buggy due to the delay and the cursor can be on other parts of the map at the time is displayed.


texts: How to calculate the impact and extent of damage of a potential tailings dam breach with ArcGIS?

texts:
I would like to be able to delineate an approximate area of a potential tailings dam breach for an emergency response plan purposes. The tailins dam and pond have not yet been built.

The data I have: LiDAR DEM (1m resolution) and provincial DEM (20m resolution), 3D drawings of the proposed tailings dam, the maximum elevation of the future tailings pond surface.

What I can do with the existing data: model the 3D drawing into the existing LiDAR surface to produce an "after" DEM/TIN surface and obtain a relatively accurate volume of the tailings at any point in time.

What I would like: simulate a breach at any random point along the dam and calculate the approximate boundary of the potential spill extent downstream.

My tools: ArcGIS Desktop Standard, Spatial Analyst, 3D Analyst, AutoCAD 2015 Map 3D (Limited experience), QGIS (Limited Experience)

I don't have access to ArcInfo functionality so some Arc Hydro functionality is not available to me.


texts: ArcGIS API for Flex create customised widget

texts:
I want to create customised widget using ArcGIS API for flex (in adobe flash builder). I surfed internet and all the samples I found of widgets use ArcGIS viewer for flex. Has anybody has the experience on creating widget using adobe flash builder? where should I start?

I am using Adobe Flex SDK 4.6 and ArcGIS API for Flex 3.6, thanks.


texts: Calculating Change in Attributes from two different years on two overlaid maps

texts:
I am looking at comparing wetlands gains and losses along a certain area. I have downloaded maps from NLCD over multiple years and am struggling on how to compare them on ArcGIS.

The two aspects that I am trying to compare are listed in the attributes tables of two datasets. I tried using an intersect function to perform the task but that datasets didn't work as a possible input.

Is there a way that I can compare these data sets to create a new layer that will show the changes between them?


texts:
I have QGIS Noosa installed on two separate PCs working on Xubuntu. Both the PCs are configured for the same QGIS repository.

However, on one of the PCs, I have the following error:

Traceback (most recent call last):   File "", line 1, in    File "/usr/lib/python3/dist-packages/qgis/__init__.py", line 80, in
    import qgis.gui   File "/usr/lib/python3/dist-packages/qgis/gui/__init__.py", line 27, in
    from qgis._gui import * RuntimeError: the sip module implements API v12.0 to v12.4 but the qgis._gui module requires API v12.5

I have already tried the following:

Doing a complete purge of the qgis software, and reinstallation.
Cleaning pip cache and reinstallation, both for Python versions 3.x and 2.x

How can this issue be rectified?


texts: Create Python Dictionary from feature classes, fields, and domains

texts:
I want to create a Python dictionary with a Feature class as the key and the fields of the feature class and domains of the feature class as separate values. I'm starting off by simply attempting to add a Feature class (as a key) and its respective feature types(as values) but I'm not sure how to build that out. here's what I am trying to do by using list comprehension but my dictionary is empty:

    Dict = {fc.name: fc.featureType.values()for fc in arcpy.ListFeatureClasses(r"gdb")}

    print Dict

Any clues on how I can start this off?


texts: How to use bold text in QGIS font marker symbol?

texts: Is it possible to use bold text with the Font marker symbol in a QGIS (point) layer?


texts: Can MakeFeatureLayer Object be Passed into ExportReport Function Layer Parameter?

texts:
I'm trying to pass an in memory layer object created from a feature class using
Code:

arcpy.MakeFeatureLayer_management()

into the
Code:

arcpy.mapping.ExportReport()

function. It's erroring out, so I think it's having trouble with the in memory feature layer (Does a feature layer have a.lyr extension?).

My code is:

Code:

arcpy.MakeFeatureLayer_management("SOILFC", "soil.lyr")
arcpy.mapping.ExportReport("soil.lyr", reportfilename, r"C:\Test_Outputs\ProjectReport1.xls")

I'm not sure if the mapping class needs an actual map document to work on, but I'm hoping not. I don't want to go through the process of creating a temp .mxd, setting a dataframe reference and then add in some layers somewhere just for the purpose of getting a layer object to create a report. I was just hoping I could create the report using in memory objects because the only thing I care about if the .xls report from the final output.


texts: Discrepancy in conversion from NAD27 to NAD83

texts:
In Southern California, where I live, multiple sources of information seem to say that the conversion between NAD27 and NAD83 is almost purely in the east-west direction (about 80 meters), with almost no correction (only a few meters) along north-south. For example, the gdal utility cs2cs seems to verify this, by showing almost no change in the latitude and a much larger change in the longitude:

$ echo "33.760268  -116.682251" | cs2cs -rsf "%.6f" +proj=latlong +datum=NAD27 +to +proj=latlong +zone=11 +datum=NAD83 -
33.760288   -116.683102 0.000032

But I'm running into weird discrepancies that I don't understand. Here are the bottom left corners of two USGS topos, one from 1988 (left) and one from 2018 (right):

The text in the 1988 topo says:

  1927 North American Datum.
  To place on the predicted North American Datum 1983 move the projection lines 80 meters east as shown by dashed corner ticks.

However, there is a clearly visible additional shift of about 200 meters in the y direction. For example, the y=35 km line coincides with a road in the old topo, but is far south of the same road in the new one.

I'm also seeing this with a handheld GPS unit. If I set it to read out in NAD27 and take the UTM coordinates on my front porch, I get this (x y):

0410910 3750091

If I then set it to NAD83 and repeat, I get this:

0410830 3750288

The sign and magnitude of the x offset are exactly as calculated by cs2cs and as claimed on the topo maps. However, there is an additional y offset of almost 200 meters.

What am I not understanding here?

If anyone can help me with this, feel free to come by my UTM coordinates for a beer. Just make sure to resolve the ambiguity first, or you'll end up at my neighbor's house up the street to the north.

[EDIT] It seems that each NAD datum has its own conversion between lat-lon and UTM:

$ echo "33.760268  -116.682251" | cs2cs -rf "%.6f" +proj=latlong +datum=NAD83 +to +proj=utm +zone=11 +datum=NAD83 -
529425.663755   3735620.830534 0.000000
$ echo "33.760268  -116.682251" | cs2cs -rf "%.6f" +proj=latlong +datum=NAD27 +to +proj=utm +zone=11 +datum=NAD27 -
529426.322142   3735427.157052 0.000000

This can also be seen on the USGS maps, where the corners of the map are defined by lat-lon, but the UTM grid is shifted north-south relative to the corners.

So I guess this sort of clears it up for me, in the sense that it shows there are really four different coordinate systems: (NAD27,lat-lon), (NAD27,utm), (NAD83,lat-lon), and (NAD83,utm). I would still be happy to get an answer explaining this further. Are there two things that have to be specified, an ellipsoid and a projection?


texts: Unable to install geopandas using conda install

texts:
I'm trying to install GeoPandas using the Command Prompt terminal. I have Python 3 and ArcGIS installed. I keep getting these messages that say there is a conflict with various packages. For example:
Package jupyterlab_launcher conflicts for:
jupyterlab -> jupyterlab_launcher[version='<0.2|>=0.10.0,<0.11.0|>=0.11.2,<0.12.0|>=0.13.1,<0.14.0|>=0.6.0,<0.7.0|>=0.5.4,<0.6.0|>=0.5.4|>=0.4.0|>=0.3.0|>=0.2.3|>=0.11.0,<0.12.0']
_ipyw_jlab_nb_ext_conf -> jupyterlab -> jupyterlab_launcher[version='<0.2|>=0.10.0,<0.11.0|>=0.11.2,<0.12.0|>=0.13.1,<0.14.0|>=0.6.0,<0.7.0|>=0.5.4,<0.6.0|>=0.5.4|>=0.4.0|>=0.3.0|>=0.2.3|>=0.11.0,<0.12.0']
I used both conda install -c conda-forge geopandas and conda install --channel conda-forge geopandas.


texts: Adding Eastings and Northings as attributes automatically to a point layer in QGIS

texts:
I have a point layer that I am adding to and am working in British National Grid. Two of the columns are Easting and Northing. When I add a new point I want these two columns to be filled out automatically with the Easting and Northing values.

Does anyone know how to achieve this using QGIS?


texts: How to define styles for features in JSON file

texts:
Styles for rendering features in OpenLayers 3 are defined as:

new ol.style.Style({
  image: new ol.style.Circle({
    radius: 6,
    stroke: new ol.style.Stroke({
      color: 'white',
      width: 2
    }),
    fill: new ol.style.Fill({
      color: 'green'
    })
  })
})

I would like to define my styles in JSON file, not in code. Is there any way how to load styles from JSON? I could easily load attributes like "color" : "green", but I cannot load something like new ol.style.Fill...


texts: TypeError: 'int' object is not iterable in reclassifying raster using arcpy

texts:
I am trying to reclassify a raster, using arcpy in Python. My code is below

from osgeo import gdal
import arcpy
from arcpy import env
from arcpy.sa import *

raster_file = 'M:/tgopalakrishna\Data\CONUS_def_001_in_020\CONUS_def_001_in_020.tif'

#Opening the raster file
dataset = gdal.Open(raster_file)
band = dataset.GetRasterBand(1)
print band

#Reading the raster properties
projectionfrom = dataset.GetProjection()
print projectionfrom
geotransform = dataset.GetGeoTransform()
print geotransform

block_sizes = band.GetBlockSize()
x_block_size = block_sizes[0]
print x_block_size
y_block_size = block_sizes[1]
print y_block_size

xsize = band.XSize
print xsize
ysize = band.YSize
print ysize

max_value = band.GetMaximum()
print max_value
min_value = band.GetMinimum()
print min_value

stats= band.GetStatistics(0,1)
print stats # prints min, max, mean and std

# Set local variables
reclassField = "VALUE"
remap = RemapRange([0,999,1])

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

outReclassify = Reclassify(raster_file, reclassField, remap, "NODATA")

outReclassify.save("M:/tgopalakrishna\Outputs\Defended_state/trial.tif")

The errors I get are

  Traceback (most recent call last):   File "*************trial.py",
  line 41, in
      remap = RemapRange([0,999,1])

  File "D:\Program Files
  (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\ParameterClasses.py", line
  1160, in init
      _Remap.init(self, remapTable)

  File "D:\Program Files
  (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\ParameterClasses.py", line
  1129, in init
      self.remapTable = [list(record) for record in remapTable]

  TypeError: 'int' object is not iterable

So clearly, there is something wrong with the range of the reclassification. How do I solve this?


texts: How to write the results of lidR grid_metrics to a georeferenced raster?

texts:
I'm using the lidR package to calculate LiDAR grid metrics (p.18 documentation) for forestry applications. For example, the following calculates the mean elevation within a 20m cell:

library(lidR)

LASfile <- "/Users/aaron/Desktop/data/1958-08-53.laz"
lidar <- readLAS(LASfile)

l <- grid_metrics(lidar, mean(Z), 20)

I can plot the results using the following approach:

grid_metrics(lidar, mean(Z), 20) %>% plot

What is the proper method to write the grid metrics to a georeferenced raster (e.g. geotiff)?


texts: Adding certain number for each community with QGIS field calculator?

texts:
I have a big data with the parcel numbers of different communities. Now  I want to join the parcel numbers of the attribute list with an excel table. Since they need to have something in common to realize the join, I have to get the parcel numbers into the same structure like the excel table.

The parcel number of the attribute list look for example like this: [803, 1]
In my excel table I created the same structure, but since there might be the same number in different communities, I had to separate them, that’s why I choose to use another number for each community.
For example: [803, 1]_1 (1 for community A) , [803, 1]_2 (2 for community B).

The problem is that I don´t know how to create this formula in my attribute table without having to edit each parcel number. Is there a formula to use? In the attribute table I have the parcel number like this [803, 1] and also the community number like this [08, 8030]. I add a picture of my attribute table, so I think its easier to understand what I need.

Is there a formula that says for example something like this: when in “Gemarkung” [08, 8030], then in “Flurstueck” [621, 4]_1 and when in “Gemarkung” [08, 8015], then in “Flurstueck” [195, ]_2 (of course in "Flurstueck should be the same number and just add _2 or _1...) I hope I could explain myself with what I need.

I am using QGIS version 2.12.3 Lyon.


texts: Basemap error offset for QuickMapServices QGIS 3.0.0. in Windows

texts: I noticed an offset error for all basemaps provided by Quickmapservices in Qgis 3.0.0 installed in Windows. It was installed with the setup.exe application (advanced options, 4 GB after installed). The offset error is of about 47 km measured from north to south. Does anyone notice the same error and have an easy fix for it? I reinstalled this plugin several times but the error persists. Any solution is very welcomed.


texts: Using QGIS 2.12 Lyon on a Mac (Intel, 2013). I check the SPIT line in plug-ins, but the plugin does not appear in the menu bar as promised. Log out, restart, nothing helps. Other plug-ins work fine and appear in the bar. Ideas?


texts: How to "output to a raster file" in a standalone app based on a composition?

texts:
In a standalone app I want to create a print composer with the map, a legend and a scale bar. I try with different codes but without success. I search in several foruns but I didn't find an answer.

 basename_R = fileInfo_R.baseName()
 R_color_raster = QgsRasterLayer(R_color, basename_R)

 layerset = []
 QgsMapLayerRegistry.instance().addMapLayer(R_color_raster)
 layerset.append(R_color_raster.id())
 myMap = QgsMapSettings()

 myComposition = QgsComposition(myMap)
 myComposition.setPlotStyle(QgsComposition.Print)

 #legend
 legend = QgsComposerLegend(myComposition)
 legend.model().setLayerSet(myMap.layers())
 legend.model().setLayerSet
 myComposition.addComposerLegend(legend)
 legend.setItemPosition (25,122)
 legend.setFrameEnabled(True)
 legend.setScale(.7)

 output_jpeg = "..../Factor.jpg"
 myImage = myComposition.printPageAsRaster(0)
 myImage.save(output_jpeg)

Legend works fine but I have no success with map composer. I try everything but nothing. Any ideas?

UPDATE

So I find this information that I think its useful for anyone that uses this composer option: Making serial maps from template using PyQGIS?

In version 2.4 and above the maprenderer is replaced by mapsettings. I try this and I am still trying but without success in compose my map. Does anyone try and get this right?

UPDATE

Ok, I find the solution to "add map". Here is the code:

fileInfo_R = QFileInfo(R_color)
 basename_R = fileInfo_R.baseName()
 path = fileInfo_R.filePath()
 R_color_raster = QgsRasterLayer(path, basename_R)

 layerset = []
 QgsMapLayerRegistry.instance().addMapLayer(R_color_raster)
 layerset.append(R_color_raster.id())
 myMap = QgsMapSettings()

 mapRectangle = myMap.fullExtent()
 myMap.setExtent(mapRectangle)
 myMap.setOutputSize(QSize(1300,900))

 myComposition = QgsComposition(myMap)
 myComposition.setPlotStyle(QgsComposition.Print)

 x,y = 0,0
 w,h = myComposition.paperWidth(), myComposition.paperHeight()
 composerMap = QgsComposerMap(myComposition,x,y,w,h)
 composerMap.gridEnabled()
 myComposition.addComposerMap(composerMap)
 composerMap.setLayerSet(myMap.layers())

 #legend
 legend = QgsComposerLegend(myComposition)
 legend.model().setLayerSet(myMap.layers())
 legend.model().setLayerSet
 myComposition.addComposerLegend(legend)
 legend.setItemPosition (25,122)
 legend.setFrameEnabled(True)
 legend.setScale(.7)

 item = QgsComposerScaleBar(myComposition)
 myComposition.addComposerScaleBar(item)
 item.setStyle('Numeric')
 item.setItemPosition(0,122)
 item.update()

 output_jpeg = "..../Factor.jpg"
 myImage = myComposition.printPageAsRaster(0)
 myImage.save(output_jpeg)

Ok, but now I don't need the legend, I need to insert a picture and I try with addComposerPicture but without success. Any idea? Also I want to add the scale bar but the numeric gives me 1:1. I assume that the coordinate system isn't readed. Any ideas?


texts:
I have an OpenLayers vector layer that loads data form service. But it loads all data initially.

        var mylayer = new OpenLayers.Layer.Vector(null, {
            "protocol": myWFSProtocol,
            "styleMap": style,
            "strategies": [
                new OpenLayers.Strategy.BBOX(),
                new OpenLayers.Strategy.Cluster()
            ]
            "filter": new OpenLayers.Filter.Logical({
                type: OpenLayers.Filter.Logical.OR,
                filters: [
                    new OpenLayers.Filter.Comparison({
                        "type": OpenLayers.Filter.Comparison.NOT_EQUAL_TO,
                        "property": attributeNames.catid,
                        "value": "null"
                    })
                ]
            })
        });

My filter says that get data if catid is not null, so so many data is loading at application initialized.

I want do not want to initialize it at first. I want to get data on button click. Is there any property to initializing Layer or protocol.


texts: SQL-Like Inner Join on IDs in Two Layers in QGIS

texts:
Using QGIS v2.6.0-Brighton.

This question is not about editing layers geographically in anyway. It is merely trying to find data from two layers as you would in an SQL Query and append an extra ID to the second layer.

I have one polyline layer each identified with two different IDs. There is a point layer also with one of these IDs, I wish to find this corresponding ID and append them as a column to the point layer.

How do I do this without finding each row individually and typing it in manually (i.e. larger scale)? I would prefer to use SQL like syntax but if there is a GUI as well that would be ideal. Happy to install official plugins.


texts: Creating lines from polygon layers boundaries?

texts:
I need to find boundaries between two polygon layers and express them as lines. There are plenty of polygons in each layer, so I'd like to avoid doing it "by hand".

Example:

Red line shows how boundaries should look.


texts: Unable to update integers in PGDB table field using ArcObjects

texts:
Dim pLinerow As IRow
pLinerow = ptable.CreateRow
pLinerow.Value(pLinerow.Fields.FindField("IS_CLOSED")) = 0

when this code is executed,  the below mentioned error is returned,

  "Error HRESULT E_FAIL has been returned from a call to a COM
  component"


texts: ArcMap point symbol/marker size, How it works?

texts:
I am trying to figure out how the size of markers in arcgis works.
I understand that a square marker set to 1 pt is 1 pt for specific scale. but when I import emf or font, How does the change size works then? if like the symbol is rectangular or a picture of a bus?

If we take the bus, I set 20 pt size does it take 20 point of the smallest part of the picture or what? Please help! :)

And also, how can I set a specific size to be correct for illustrator export? the size gets different, is there any script for arcmap or illustrator to set correct size of symbols?


texts: Arcpy.da.SearchCursor() help when attribute doesn't exist

texts:
I'm updating text elements off of an arcpy.da.SearchCursor script and needed to construct the script to return a response for when the attribute does NOT exist. Here is the snippet of the code THAT WORKS:

fcWO = r"W:\\GIS_Projects\\Impervious_Surfaces_Tables\\SUM_WORKORDER"
field1 = "SURFACE_TYPE"
field2 = "SUM_SQ_FT"

valueB = r"N/A"
valueM = r"N/A"
valueP = r"N/A"
valueS = r"N/A"
valueR = r"N/A"

with arcpy.da.SearchCursor(fcWO, (field1, field2)) as cursor:
    for row in cursor:
        if row[0] == "BUILDING":
            valueB = row[1]
        elif row[0] == "MISC":
            valueM = row[1]
        elif row[0] == "PARKING":
            valueP = row[1]
        elif row[0] == "ROAD":
            valueR = row[1]
        elif row[0] == "SIDEWALK":
            valueS = row[1]

for elm in arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT"):
    if elm.name == 'Building':
        elm.text = valueB
    elif elm.name == 'Misc':
        elm.text = valueM
    elif elm.name == 'Parking':
        elm.text = valueP
    elif elm.name == 'Roadway':
        elm.text = valueR
    elif elm.name == 'Sidwalk':
        elm.text = valueS


texts: When generate buffer for zonal statistic, how to avoid the neighbor object?

texts:
I have a layer of many polygons indicating water area, scattered everywhere. Some of them are very close to another.

I want to know what's the slop around the water sites. So I already have a layer of slope all across the whole area.

I intended to generate 50m, 100, 150 buffer, then use the tool "Zonal statistic as table" to give me the slope around those waters' data. But I found some polygons are located in the buffer area of another polygon. As you may guess, slope data on the water site is 0. If I go to zonal analysis directly I will get a wrong answer about the average slope.

So I want to eliminate the water area located in the buffer zone.

How can I do that using ArcGIS Desktop 10.3?


texts:
Does maptiler support DEM *.tif images?

I was generating terrain tiles via cesium-terrain-builder but I wonder can I use maptiler somehow?


texts: How can I generate PDF file in python on QGIS?

texts:
I am following these instructions:

https://bitbucket.org/rptlab/reportlab

but i am getting this error:

Couldn't load plugin Importtool due to an error when calling its classFactory() method

Traceback (most recent call last):

  File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 219, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/n/.qgis2/python/plugins\Importtool\__init__.py", line 34, in classFactory
    from .Import_tool import Importtool

  File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 478, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)

  File "C:/Users/n/.qgis2/python/plugins\Importtool\Import_tool.py", line 39, in
    import reportlab

  File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 478, in _import

    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named reportlab


texts: I made a color composite using Raster > Miscellaneous > Build Virtual Raster in QGIS. The output is a VRT file but I need TIFF. Then I went through Raster > Conversion > Translate. The output is TIFF however the software automatically change the order of bands and then the output TIFF file has different color. How can I fix that?


texts: Leaflet opacity slider with L.tileLayer.wms

texts:
I would like to control the opacity of a wms layer (using L.tileLayer.wms) in Leaflet using an opacity slider.
In L.tileLayer.wms I can control the opacity by setting the opacity:value.
This is my piece of code; The wms layer
var wmsLayer_2 = L.tileLayer.wms('https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php', {
  layers: ['CP.CadastralParcel'],
  crs: crs_6706,
    format: 'image/png',
  maxZoom: 21,
  transparent: true,
}).addTo(map);

The slider implemented in the leaflet control
function updateOpacity(value) {
        wmsLayer_2.setParams({opacity:value});
        };

var control = L.control.layers(baseMaps,
    {'Opacity_layer <br><input id="range" type="range" min="0" max="1" step="0.1" value="1" class="slider" onchange="updateOpacity(this.value)">':wmsLayer_2},{collapsed: false}).addTo(map);

However I can't manage the opacity.
Here a complete example


texts: Deleting datastore (and its content) with GeoServer REST API?

texts:
I need to delete a Datastore (as well as all its content) using GeoServer's REST API.

I'm trying:

curl -u usr:pwd -XDELETE http://example.com/geoserver/rest/workspaces/ws/datastores/my_datastore.xml

But the response I get is:

  datastore not empty

So, I assume I need to delete all layers in My Datastore. Fine.

Deleting all layers in the datastore:

I couldn't find a way to list only the layers in a concrete DS, so I did:

curl -u usr:pwd -XGET http://example.com/geoserver/rest/layers.json

that returns something like:

{
    "layers":{
        "layer":[
            ..., // Some unrelated layers
            {
                "name":"my_datastore.poi",
                "href":"http://example.com/geoserver/layers/my_datastore.poi.json"
            },
            {
                "name":"my_datastore.sector",
                "href":"http://example.com/geoserver/layers/my_datastore.sector.json"
            },
            ... // More unrelated layers
        ]
    }
}

After filtering this response to get only the layers in the datastore to be deleted, I delete every layer doing:

curl -u usr:pwd -XDELETE http://example.com/geoserver/layers/my_datastore.poi.json
curl -u usr:pwd -XDELETE http://example.com/geoserver/layers/my_datastore.sector.json

And it apparently works (if I check in the GeoServer's web interface, the layers are gone...)

Finally, deleting the supposedly empty datastore

After deleting the layers I try again:

curl -u usr:pwd -XDELETE http://example.com/geoserver/rest/workspaces/ws/datastores/my_datastore.xml

But:

  datastore not empty

So, my double question is:

What else do I need to delete for the datastore to be empty?
Is there any easier way to delete a Datastore (and its content) via REST API?


texts: Arcpy ListFeatureClasses in same workspace differs when run from ArcMap VS IDLE

texts:
This Works

The following python script works when connecting to my dev ArcGIS server, from both the ArcMap python window and the IDLE python GUI/shell running on my dev workstation:

import arcpy
arcpy.env.workspace = "Database Connections/Connection to 172.20.0.XX.sde"
str(arcpy.ListFeatureClasses("*"))

Works: Prints list of all the feature classes in the Enterprise GeoDB

This Does NOT Work

The following python script, which is the same as above but aside from the env.workspace value, DOES NOT work when connecting to my prod ArcGIS server from the IDLE python GUI/shell running on that same ArcGIS server, but DOES work from the ArcMap python window on that server...

When run from ArcMap 10.2.1 python window:

import arcpy
arcpy.env.workspace = "Database Connections/Connection to 172.20.0.XX.sde"
str(arcpy.ListFeatureClasses("*"))

Works: Prints list of all the feature classes in the Enterprise GeoDB

When run from IDLE python GUI/shell:

import arcpy
arcpy.env.workspace = "Database Connections/Connection to 172.20.0.XX.sde"
str(arcpy.ListFeatureClasses("*"))

Fails: Prints 'None'

The Windows user I'm logged into the prod ArcGIS server with is an admin, but The only thing I can come up with is

Notes

Note I have tried all of the following variations when setting the workspace:

arcpy.env.workspace = r"Database Connections\Connection to 172.20.0.49.sde"
arcpy.env.workspace = "Database Connections\Connection to 172.20.0.49.sde"
arcpy.env.workspace = "Database Connections" + os.sep + "Connection to 172.20.0.49.sde"
arcpy.env.workspace = os.path.join("Database Connections", "Connection to 172.20.0.49.sde")
Etc...


texts:
After training the kmeans model, how can i obtain the number of classes of the model? i've tried with get(), getInfo() and doesnt work.

https://code.earthengine.google.com/ff0c7bd8102d17b2d267b96bd070342a

Value that must to be extracted

// STEP 1: Import image we wish to classify
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1')

// Remove the clouds
var cloud_free = ee.Algorithms.Landsat.simpleComposite({
  collection: l8.filterDate('2018-01-01', '2018-12-31'),
  asFloat: true
})

// STEP 2: Draw your box which you will sample from!
// Make the training dataset.
var training = cloud_free.sample({
  region: geometry,
  scale: 30,
  numPixels: 5000
});

// Step 3: Create the clusterer and train it, play around with
// number of classes
var clusterer = ee.Clusterer.wekaKMeans(3).train(training);

print('clusterer',clusterer instanceof ee.ComputedObject);

var to_print = clusterer.getInfo();

print('clusterer_info',clusterer.toString());

print(typeof clusterer);

var clusterer_list = ee.List(clusterer);

print('clusterer_list',clusterer_list);

print(typeof clusterer_list);

//var to_print_2 = clusterer.get('nClusters');

//print(to_print_2);

//var to_print_3 = clusterer.Filter.filterMetadata('nClusters','equals',3);

//print(to_print_3);

// Step 4: Classify our image
var classified = cloud_free.cluster(clusterer);
print(classified);
var cluster_1 = classified.eq(0);
print('min',cluster_1);

// Display
Map.setCenter(-122.31040963183364, 37.782065864682096, 11)
Map.addLayer(classified.randomVisualizer(), {}, 'clusters');


texts: Arc 10.2.2 crashes when adding query layer?

texts:
I'm trying to add a query layer from an Oracle database into a blank 10.2.2. map.   I enter the query, validate, finish, and the software thinks for a second and then just winks out.  No error message, no nothing.

One of my co-workers can perform the exact same query layer and a table shows up in his blank map.  The only difference is that we are working in different paths on the network drive and we have some slightly different software installed on our computers,  (I have InfoMaker, for example).  I am able to perform the queries from InfoMaker and Oracle Developer, just not ArcGIS.

We both have the same model of computer (Dell 7010) and the same Zen image installation of the software.  I deleted normal.mxt.


texts: Export to drive - User memory limit exceeded

texts:
I am using the GSW seasonality raster in GEE. I started by clipping the GSW to the whole South America region and masked this region only. However, I do not manage to export to GDrive and I get the error:

  User memory limit exceeded

Here is the code:

//South America polygon
SA_poly = ee.FeatureCollection("South_America_polygon");

//GSW 1.1 with seasonality and clip to South America shoreline
var seasonality = ee.Image("JRC/GSW1_1/GlobalSurfaceWater")
                   .select('seasonality').clip(SA_poly);

//export to drive
Export.image.toDrive({
  image: seasonality,
  description: 'seasonality',
  scale: 30,
 folder: 'GEE',
  fileFormat: 'GeoTIFF',
  region: SA_poly,
  crs: 'EPSG:4326',
  skipEmptyTiles: true,
});

Is there a way to overcome that? Create a split polygon function of the export polygon to a size that would allow GEE export?


texts: Building footprints of the city of Guayaquil, Ecuador

texts:
Does anyone have/know where I can find a shapefile/dataset of the building footprints of the city of Guayaquil in Ecuador?

this type of information is so hard to get here in Ecuador. I went through all of the government database pages but there is none containing the information i am looking for.


texts:
Using the tz_world.shp file, I want to find the polygon in which a point falls.  Technical details - language: R, v3.2.3; platform: Windows 7 professional (64-bit, quad core, 32 GB RAM).

I have 12169 latitude/longitude pairs and many polygons to search.  Using the bounding box for each polygon against each entry in the 12169 x 2 double array of longitude/latitude points would reduce the search time.  I have so far:

timeZonesShpFile = readShapeSpatial( 'tz_world.shp' )
timeZoneList = split( timeZonesShpFile, timeZonesShpFile$TZID )
ps = lapply( timeZoneList, Polygon )
p1 = lapply(seq_along(ps), function(i) Polygons(list(ps[[i]]), ID = names( timeZoneList )[i] ) )
my_spatial_polys = SpatialPolygons( p1, proj4string = CRS("+proj=longlat +datum=WGS84") )

>my_spatial_polys[1]@bbox
    min       max
x -6.091949 -5.553957
y  4.980661  7.623367

Using extract() is very slow: value = extract( timeZonesShpFile, data.frame( lon = aptLon[i], lat = aptLat[i] ) )

I think that the fast solution is to reduce the dimensionality of the search by keeping only those bounding boxes (and associated polygon indices) in which the ith latitude/longitude point value falls.  Also, occasional records in the latitude/longitude array have the values of "#VALUE!", which I can catch with is.numeric(unlist())==TRUE.

Given the bounding box, I can test latPoint >= min(latBBox) & latPoint <= max(latBBox) & lonPoint >= min(lonBBox) & lonPoint <= max(lonBBox).  I can create a shorter listing of polygons against which to use extract() for the final test.

In summary, I want to inspect the bounding box for each polygon, retain the polygon index, and test each latitude/longitude location within each box, thereby retaining the list of polygons in which the ith location falls.


texts: Metric/parameters to compare length of polyline

texts:
What would be a good metric to compare the length of two polylines?

A very naive approach would be to define a maximum delta which the two lines can differ.

I'm looking for more approaches. For example one that might also take a function over the length of the polyline as a parameter. Does anybody has algorithms/functions to share and has even evaluated parameters that work quite good?


texts: Rotating a selection of features around their individual centroid rather than the combined centroid

texts:
The rotate feature in QGIS 3 allows you to select multiple features at once and rotate them around the centroid of their combined geometries, as found in the documentation.
You can also use an anchor point different from the default feature centroid by pressing Ctrl or if you hold Shift before clicking on the map, the rotation will be done in 45 degree steps, which can be modified afterwards in the user input widget.
I would like to rotate a selection of multiple features around each features individual centroid?
Is this possible?


texts: Deleting feature classes with similar name in ArcPy?

texts:
I have a script that does a dissolve on polygons with a common field. I'm only interested in keeping the outputs that contain more than one polygon. I'm new to python and having trouble writing a script that will delete all the polygons that have a name starting with 'dissolve_*'. I get an error saying "name 'fnmatch' is not defined. Below is the script. The problem is at the very end of it.

Found the Solution  I had to move the Delete_management outside of the SearchCursor. I've updated the scritp to one that works and deleted extra bits as per PolyGeos request.

I'm using ArcGIS 10.2 Desktop

import arcpy, os, sys
arcpy.env.overwriteOutput = True
arcpy.env.workspace = arcpy.GetParameterAsText(0)
PD_ED = arcpy.GetParameterAsText(1)
NullwhereClause = arcpy.GetParameterAsText(2)
output_dir = arcpy.GetParameterAsText(3)

# Local variables for dissolve:
PD_ED_fld = ['Add_PNTID', 'PD_NO', 'ED_NO']
PD_PNTID_Lay = 'PD_PNTID_Lay'
dissolveField = ['Add_PNTID']
dissolve_lay ="lay_dis"
try:
    with arcpy.da.SearchCursor(PD_ED, PD_ED_fld, TestPDwhereClause) as PDcursor:
        for PDrow in PDcursor:
            # Returns string values for the field names specified
            PDPNTIDrow = str(PDrow[0])

            # where statement to select a single PNTID value from the field
            PDwhereClause = 'Add_PNTID =' +PDPNTIDrow

            # Process: Make Feature Layer
            arcpy.MakeFeatureLayer_management(PD_ED, PD_PNTID_Lay, PDwhereClause)

            dissolve = arcpy.Dissolve_management(PD_PNTID_Lay, disPDoutput, dissolveField, "", "SINGLE_PART", "")

            # Process: Make Feature Layer
            dissolveMakLay = arcpy.MakeFeatureLayer_management(dissolve, dissolve_lay)
            # Process: Get Count - if the count in the layer is greater than 1 an output is created otherwise "Only one polygon found" is printed
            dissolveCount = int(arcpy.GetCount_management(dissolveMakLay).getOutput(0))
            if dissolveCount > 1:
                arcpy.FeatureClassToFeatureClass_conversion(dissolveMakLay, output_dir, disPDoutManyPoly)
                print("More than one Polygon found")
            else:
                print("Only one polygon found")

except arcpy.ExecuteError:
    print arcpy.GetMessages(2)
    arcpy.AddError(arcpy.GetMessages(2))
except Exception as e:
    print e.args[0]
    arcpy.AddError(e.args[0])

try:
    for fc in arcpy.ListFeatureClasses("dissolve_*"):
        arcpy.Delete_management(fc)

except arcpy.ExecuteError:
    print arcpy.GetMessages(2)
    arcpy.AddError(arcpy.GetMessages(2))
except Exception as e:
    print e.args[0]
    arcpy.AddError(e.args[0])


texts: Using pgrouting and algorithms in PostgreSQL?

texts:
I am trying my first pgrouting in a roads shapefile I load it in PostgreSQL first, the table has osm_id, gid, name, ref, type, oneway, maxspeed

Then I add columns source and target in table:

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

Then I make the topology with function:

SELECT pgr_createTopology('roads', 0.0001, 'geom', 'gid');

make indexis

CREATE INDEX roads_source_idx ON roads("source");
CREATE INDEX roads_target_idx ON roads("target");

Add cost column and count the cost with geom column use ST_LENGTH method

ALTER TABLE roads ADD COLUMN length double precision;
UPDATE roads SET length = ST_LENGTH(geom);

Add reverse cost column and count it

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

Finally I run the distrack query, which is running in PostgreSQL, but if I try to run this query in QGIS it is not running. This is the first problem.

SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
                SELECT gid AS id,
                         source::integer,
                         target::integer,
                         shape_leng::double precision AS cost
                        FROM roads',
                30, 300, false, false);

Then I try to use A*

I count the x1 y1 x2 y2 but the columns is empty and I don't know why
insert the columns with these queries:

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

Then I use these queries to field them but the columns still empty after execution:

UPDATE roads SET x1 = ST_x(ST_startpoint(geom));
UPDATE roads SET y1 = ST_y(ST_startpoint(geom));
UPDATE roads SET x2 = ST_x(ST_endpoint(geom));
UPDATE roads SET y2 = ST_y(ST_endpoint(geom));

I read in net same versions have problem with startpoint and endpoint and I use the PointN but I still take empty columns:

UPDATE roads SET x1 = ST_x(ST_PointN(geom, 1));
UPDATE roads SET y1 = ST_y(ST_PointN(geom, 1));
UPDATE roads SET x2 = ST_x(ST_PointN(geom, ST_NumPoints(geom)));
UPDATE roads SET y2 = ST_y(ST_PointN(geom, ST_NumPoints(geom)));

And the query for A*

SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_astar('
                SELECT gid AS id,
                         source::integer,
                         target::integer,
                         shape_leng::double precision AS cost,
                         x1, y1, x2, y2
                        FROM roads',
                30, 60, false, false);

How can I show the routing in QGIS from PostgreSQL I did it, but I still
have problem with A* x1,y1 & x2,y2 after the update is still empty the  columns x1 y1 x2 y2?

I try to use these queries but the x2y2 is not filling:

UPDATE roads SET x1 = ST_X(ST_StartPoint(ST_geometryN(geom,1)));
UPDATE roads SET y1 = ST_Y(ST_StartPoint(ST_geometryN(geom,1)));
UPDATE roads SET x2 = ST_X(ST_EndPoint(ST_geometryN(geom,ST_NumPoints(geom))));
UPDATE roads SET y2 = ST_Y(ST_EndPoint(ST_geometryN(geom,ST_NumPoints(geom))));

Please if anyone can help me


texts: Loop though records in PostGIS (using PL/PYTHON or PGSQL)

texts:
I've seen many example of loops, but I can't get that one to work.

I work for a regional transit agency and each year, we conduct several surveys. I geocode everything to the house level for internal use, but I've been asked to anonymize our results so I could be sent to our partners while preserving the anonymity of the persons who answered.

We have an address point table containing every street addresses on our territory. For each of them, we have the number of housing units.

To do it, I use PostGIS. For each geocoded address, I select the 50 nearest points, and then select base on the weight of the number of units. Basically, it's weighted random n nearest neighbor selection.

If I do my script for a single point, it's perfect. Maybe it's not the most elegant solution, but it works. I've made several tests and the weighted random select from the selected 50 nearest points is good.

The only thing is that I can't do this using a function (didn't say it can't be done, just that I can't ;) )

So, I created a python script, that loops using a serial value field and it's working. The problem is that I'm the only one that uses python so I wanted to do it entirely within PostGIS.

Here's my script:

conn = psycopg2.connect("dbname=anon user=postgres password=######")

#-- Cursor
cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)

#-- Add necessary fields
cur.execute("""alter table od2008 add column val serial;
               alter table od2008 add column aqidadr varchar(20);
               alter table od2008 add column xanon numeric;
               alter table od2008 add column yanon numeric;
               alter table od2008 add column distance numeric;
               alter table od2008 add column sm2 numeric;""")

#-- Number of rows
cur.execute("SELECT max(val) as maxval FROM od2008")
rec = cur.fetchone()
nb = rec['maxval']

val = 1
while val <= nb:

    #-- Select row "i"
    cur.execute("SELECT * FROM od2008 WHERE val = %i" %(val))
    rec = cur.fetchone()

    #-- (1) Select the 50 nearest points
    #-- (2) Calculating total weight
    #-- (3) Weighted random selection

    #-- Table od2008 --> to anonymize
    #- Table aq_adr --> Address points to select from

    cur.execute("""
       -- (1)
       CREATE TABLE near50 AS
        SELECT aq_adr.*, ST_Distance(od2008.geom, aq_adr.geom) AS distance
        FROM aq_adr, od2008
            WHERE ST_DWithin(od2008.geom, aq_adr.geom, 1000) AND
            od2008.val = %i
            AND od2008.sm = aq_adr.sm ORDER BY distance LIMIT 50;

       -- (2)
   CREATE TABLE weight AS
            SELECT sum(nbunit) AS weight
        FROM near50;

       --(3)
       UPDATE od2008
            SET aqidadr = (SELECT idadr
                FROM near50, weight
                ORDER BY Random()*(1/(nbunit*100/weight.weight)) LIMIT 1)
                WHERE od2008.val = %i;
   """%(val, val))
    val += 1

    #-- Delete temp tables
    cur.execute("DROP TABLE near50, weight;")

    conn.commit()

Any tips on how I could loop through each geometry and create my 50 nearest points selection?

Thanks!!


texts: gdal_merge.py aligns incorrectly at the seam of two GeoTiffs

texts:
I am trying to merge four GeoTiffs from the National Elevation Dataset (EPSG: 4269) into a 2x2 DEM using gdal_merge.py.  For some reason, when it merges the upper two NED rasters, gdal produces a seam between the two of approx. 1 px that it fills with incorrect values (either 0 or some large negative number).  It seems like this is probably an alignment issue.  How can I eliminate this artifact?

Thanks for the help.


texts: Extracting raster information within/inside polygon (for each contour)

texts:
I have raster data of NDVI and DEM
I have already created contour from DEM. Now I'm trying to find following:

for each contour find the raster information within the contour.

Can anyone suggest available methods in GRASS?

I used r.mask and then v.what.rast it works. However, I got the ndvi values of the area that has been masked (NDVI values that lie inside the polygon plus the outside area) [image -1] https://i.stack.imgur.com/NMNlm.jpg

Actually, I need NDVI values that lie inside the contour/polygon [image-2]  https://i.stack.imgur.com/A2YWl.jpg

I may actually need to intersect each contour with the NDVI and extract the NDVI values (inside the contour).


texts:
for someone skilled, I'm sure this will be piece of cake. I'm starting implementing WMTS into own application and even that I already use WMS for years, I have some troubles with understanding WMTS.

WMTS service: http://gis.lmi.is/arcgis/rest/services/kort_isn93/kort_isn93/MapServer (Island freely available maps)
WMTS configuration: http://gis.lmi.is/arcgis/rest/services/kort_isn93/kort_isn93/MapServer/WMTS/1.0.0/WMTSCapabilities.xml

My confusions:

1. more confusion then problem

<ows:BoundingBox crs="urn:ogc:def:crs:EPSG::3057">
  <ows:LowerCorner>-22417.9302 166013.8265</ows:LowerCorner>
  <ows:UpperCorner>1048334.75710001 852747.7883</ows:UpperCorner>
</ows:BoundingBox>
<ows:WGS84BoundingBox crs="urn:ogc:def:crs:OGC:2:84">
  <ows:LowerCorner>-31.399192003034788 61.60136743056403</ows:LowerCorner>
  <ows:UpperCorner>-6.002491614773447 68.16183386054668</ows:UpperCorner>
</ows:WGS84BoundingBox>

Coordinates for lowerCorner converted from 3057 > 4326 do not match .. -28.8782279906;61.6383760341. Weird?

2. bigger problem

Every tileMatrix has defined TopLeftCorner

<TopLeftCorner>-3.45415E7 3.85236E7</TopLeftCorner>

But converted again to EPSG:4326 I get 30.6506275545;-69.2326947 . They look like just reversed, but no change in coordinates cause correct result.

For conversion I used nice web page: http://cs2cs.mygeodata.eu/ , same result anyway directly in own application.

Any hint, link to web page, whatever, is very welcome. Mainly second problem is very "problematic" as it do not allow to me compute correct bbox for map TileMatrix.


texts: Cannot load postgis query as layer in QGIS using create table command

texts:
I've a multilinestring object and a points object , where the points are located on the line. The user can set a start and end point on the line and my goal is to return a subset of the points that are within start_point and end_point on the line.

I first create a temp table (subline_table) using ST_LineSubString to get the linestring within start_point and end_point and then query my points object with ST_DWithin using this subline_table.
If I run this query in QGIS DB Manager the correct table/points are shown, but this layer will not show up as query layer once loaded.

If I store the point selection as a table, than the output is saved as a table and can be added succesfully as a postgis layer in QGIS.

My query is as follow:

-- SELECT SUBLINE WHERE INPUT POINTS ARE SNAPPED TO LINE AND STORE IN TEMP TABLE
CREATE TEMP TABLE subline_table (geom) ON COMMIT DROP AS
SELECT ST_LineSubString(line.sub,
ST_LineLocatePoint(line.sub, (SELECT geom FROM tbl_start_end AS point_start WHERE Id = 0)),
ST_LineLocatePoint(line.sub, (SELECT geom FROM tbl_start_end AS point_end WHERE Id = 1))
)
FROM (SELECT ST_LineMerge((SELECT geom FROM linestring WHERE id = 2)) AS sub) AS line;

-- SELECT POINTS WITHIN DISTANCE OF SUBLINE
-- CREATE TABLE subpoints_table AS
SELECT
  points.geom,
  points.gid,
  points.id
FROM tbl_all_points points, subline_table line
WHERE ST_DWithin(points.geom, line.geom, 20000)

Questions:

Is it possible to add a querylayer to QGIS if there is a CREATE TABLE command in the query?
If no, any suggestions how to do this query without the need to write the subline to a table?


texts: Create raster with grid cells values based on area covered by polygons

texts:
I have a polygon data with area covered by the forests (data is here - https://www.dropbox.com/s/zgckliydalljw6a/sp_data.zip?dl=0).  I want to convert polygons to raster. The value of each grid cell should be based on the area covered by polygon. For example, if grid cell size is 100m x 100m (10000m2) and polygons covered 5000m2 of this cell - the value should be 50 (or 0.5). Below, you can find my proposed solution:

Raw data:

library('rgdal')
library('raster')
library('rgeos')
library('maptools')

sp_data <- readOGR(".", "sp_data")
spplot(sp_data)

I've created new raster with 100 meters resolution and extent equal to the polygon bounding box:

r <- raster()
bb <- extent(290000, 300000, 500000, 510000)
extent(r) <- bb
res(r) <- 100

I've also created background polygon - with the values for the areas without forest equal to 0. Afterwards, I've merged these two polygon data:

sp_back <- as(extent(r), "SpatialPolygons")
proj4string(sp_back) <- proj4string(sp_data)
sp_back <- gDifference(sp_back, sp_data)
plot(sp_back, col='red')
plot(sp_data, col='green')
sp_back <- SpatialPolygonsDataFrame(sp_back, data=data.frame(value=rep(0, length(sp_back)), row.names=row.names(sp_back)))
sp_back <- spChFIDs(sp_back, "new_id")
sp_bind <- spRbind(sp_data, sp_back)
spplot(sp_bind)

Additionally, I've created new raster with reduced resolution (500 meters):

r2 <- r
res(r2) <- 500

I've rasterized polygon data to the first raster:

sp_raster <- rasterize(sp_bind, r, field="value", fun=max)
spplot(sp_raster, aspect='iso')

At the end, I've resampled the first raster values into the second raster:

sp_raster2 <- resample(sp_raster, r2, method='bilinear')
spplot(sp_raster2, aspect='iso')

Based on the results, I've have some questions:

First of all - is my solution and the result even correct?
Conversion between vector and raster somehow simplify my data (for example, loss of small polygons) and probably add some errors. The result of the resampling is affected by this simplification. Are there any alternative solution to calculate the share of polygon area in the grid cell?

Answers using R, Grass or SAGA will be welcome.


texts:
I am trying to get a snapshot of a map with Mapwingis 4.8.8.

My problem is that I tried dozen of ways but my picturebox stays desperately empty!


texts: Updating Z value within XYZ.txt file using ArcPy?

texts:
Following Creating polygon from 2 txt files (Parent ID + Coordinates) from a few weeks ago, the script works very well. Now I need to reverse the process in order to update the Z value in the original text file describing a mesh.

From a first table of Points coordinates (Nnodes, resulting from Export Feature Attribute to ASCII):

 #Nnodes
 # X Y ID Z
 4403607.48  5788195.96  56362  87.6
 4403550.41  5788192.5  56951  93.4

I want to update the Z value from Nnodes to Mesh_orig,using the Point IDs as common keys

 # Mesh_orig
 #Geom ID X Y Z
 ND 56362 4403607.48 5788195.96 86.1
 ND 56951 4403550.41 5788192.5 91.3

It prints only the last row of Nnodes instead of the whole updated mesh_dict.

 with open (r'text2.txt','r+') as Mesh:
      for row in Mesh:
          if row.startswith('ND'):
              Inline= row.split()
              mesh_dict={Inline[1]:Inline[4]}
              KeyM = mesh_dict.keys()
              ValueM = mesh_dict.values()

 with open (r'Text1.txt','r') as Npoint:
     for KeyM in mesh_dict:
         for lines in Npoint:
             Outlines= lines.split()
             point_dict={Outlines[-2]:Outlines[-1]}
             KeyP = point_dict.keys()
             ValueP = point_dict.values()
             if KeyM in point_dict:
                    point_dict.update(mesh_dict)
                    print str(KeyM)+", "+str(ValueP)

Where is it going wrong?


texts: Create Random Points for multiple features with varying points needed

texts:
I am looking for a way to randomly create points features within all National Park forest boundaries. Each boundary has a certain number of points

Example below

PARK    Count
ABLI    60    (I would need 60 randomly generated points for this)
ACAD    413   (413 randomly generated points for this...etc)
ADAM    12
AFBG    2
AGFO    62
ALFL    19

I'd like for all of this to end up in a single shapefile

I know I can use Create Random Points to do this for a single feature (aka 60 points for ABLI), but is there a good way to run this repeatedly for each park having a different count of points need?? It's over 300 features that I need points generated for, so an idea for a script, model builder, or batch code would be fantastic.

Ultimately each new point will have a unique identifier number.


texts: QGIS 2.18.1 Predominant line attribute within polygon

texts:
I've got a polygon shape (Catchments) and a lines shape (Rivers)

Each river body has an ID. What I need is to assign the ID of the predominant river within a polygon. For example, in 1 polygon there might be many rivers with different IDs ("MSPF_EM_CD" is the ID of the rivers):

As can be seen in the picture, in the attribute table of the rivers shape there is also a column with the length of each river line. What I need is to identify the longest line within each polygon, and add the ID of this longest line to the polygon shape, so I can get a new column in the polygon attribute table with the ID of the longest river within each polygon.

I am using QGIS 2.18.1 Las Palmas, in Windows 10.


texts:
Let me describe the data setup. A line has been provided as an array of coordinate pairs, one for each vertex. This array is variable in size, anywhere from 2 to 200 coordinate pairs per line. What I am trying to wrap my mind around is the best way to create these lines. Arcpy and ArcObjects are both viable options, but ArcObjects is preferred due to its speed (generally).

I have brainstormed a few methods, but they fall short for some reason or another:

1) I could use the Polyline class with arcpy to create polylines and then Copy Features, but the resulting geometry would not have any of the attribution necessary to identify said points.

2) I could use XY To Line, but this only works with one pair of point, and I require something much more robust.

The long and short is that I require a method that can create a polyline from arrays and maintain a schema for data that I have created in a table. One idea I have is to create said table and schema, and then somehow create the polyline and copy its shape property to the output table, porting over the ID and whatever other attribution is needed. How can I go about creating the polyline's shape property without resorting to an intermediate shapefile? Or, if there is a better method that comes to mind, please share it with me.

Update This code seems to work!

IFeatureClass fClass = gpUtil.OpenFeatureClassFromString(path);
IFeatureCursor cursor = fClass.Insert(true);
IFeatureBuffer buffer = fClass.CreateFeatureBuffer();

<snip...>

string[] nodearray = a.Split(';');
foreach (string i in nodearray)
{
    query.WhereClause = "\"PNT_ID\" = '" + i+"'";
    Console.WriteLine(query.WhereClause.ToString());
    ICursor search = table.Search(query, false);
    IRow row = search.NextRow();
    float x = float.Parse(row.get_Value(2).ToString());
    float y = float.Parse(row.get_Value(3).ToString());
    int id = int.Parse(row.get_Value(1).ToString());
    point.X = x;
    point.Y = y;
    point.ID = id;
    point.SpatialReference = sRef;
    pCollection.AddPoint(point, ref Missing, ref Missing);

}
IGeometry5 polygeo = (IGeometry5)pCollection;
buffer.set_Value(1, polygeo);
buffer.set_Value(2, wayid);
cursor.InsertFeature(buffer);


texts: Create shapefile based on raster values in QGIS/GDAL

texts: I need to find all pixels with a certain value and create a shapefile from those locations. Any idea how to get started ?


texts: How to Get QGIS to calculate Area in Layer's Units?

texts:
I have a Polygon Shapefile in a custom CRS which has feet as units. The Qgis Map is also set to the same CRS.

When I open the field Calculator, and enter $area, I see that the proper area is calculated:

But when I save it in the Attribute table, some other value gets saved:

The ratio between the values seems to suggest that it is doing some conversion from sqFeet to SqMeters.

I was under the impression that the Field calculator would calculate areas in the units of Layer, but here it seems to be doing something else.

Instead of manually converting, is there a way to get QGIS to calculate areas in Layer's Units?


texts: Modifying vertices x,y values to remove all value from 3rd decimal in ArcMap?

texts:
I want to check whether there is a way to update the vertices x, y value and remove all values from the 3rd decimal. I only want 2 decimals (no matter the 3rd decimal is 0 or 9, just remove it).

The original data and the desired data is like the following:

25642.882  38471.416 (need to change to 25642.88  38471.41)

25648.078  38446.104 (need to change to 25648.07  38446.10)

25655.240  38406.615 (need to change to 25655.24  38406.61)

Any advice?


texts: How to "Apply Edits" to an ArcGIS Server edit session using Python

texts:
ArcGIS Server allows the user to apply edits in an editable layer.

Eg the sample layer at http://sampleserver6.arcgisonline.com/arcgis/rest/services/Military/FeatureServer/3/applyEdits will accept a JSON input response such as [{'attributes': {'distance':1234}}]

Assuming that I have a JSON object in Python, how can I send this object to ArcGIS Server? ie, what is the Python syntax for Apply Edits?

I presume it's something along the lines of urllib2.urlopen(<endpoint>, <JSON>) but I can't figure out the exact syntax.


texts: Convert full colour height map to greyscale

texts:
I try to convert full colour map with colour range like this:

to greyscale. I use mostly Corel Photo-Paint (similiar functions to Adobe Photoshop), but maybe QGIS or ArcMap could help?


texts: which projection for preserving distances in USA (all 50 states)?

texts:
I have a shapefile of different administrative regions that are in all 50 states in the USA. I also have many lat/long coordinates, and I'm using NNJoin in QGIS to find the distance from each point to the nearest polygon.

I'm wondering which projection (if any) should I use to minimize distortion when calculating distances in the entire USA region? Or would it be better to do three batches, one for contiguous USA, one for Alaska and one for Hawaii (if so, which projections then)?


texts: How to convert a regular point grid to raster?

texts:
Possible Duplicate:
  How to convert a point grid to a raster?

I´d like to convert a regular grid (25 meter) of points (XYZ) to a raster. I have tried the GRASS-tool r.in.xyz but I always get a integer raster where all cells have a value of 1.


texts:
I'm trying to upgrade to the 10.3 version of ArgGIS (currently, I have Advanced Desktop 10.2). When I go to the download page in my ESRI account, I'm given the option to download ArcGIS for Desktop or ArcGIS Pro

Which one should I download? What's the difference between these two?

I use ArcGIS Advanced Desktop 10.2 to do the following:

network analysis using the Network Analyst tool
spatial statistics using Spatial Analyst
Lots of tasks processing raster data
Extensive python scripting using the arcpy module in python 2.7

Basically I want the one that is the fastest and can do the most while also allowing me to do all the things I currently do in 10.2


texts: Build and export PDF maps using ArcObjects?

texts:
I have a series of featureclasses and rasters generated using ArcObjects and calls to the Geoprocessor. I have the paths to the output, which are located in a file geodatabase. How can I use ArcObjects to build a map using a combination of these outputs and then save the final map as a PDF (or JPEG/PNG if PDF is not possible)? If this is possible, can I also add a legend/north arrow, etc?

The purpose of these maps is quick quality control and the output are all in the same projection and cover the same area.

I'm unfamiliar with building visual outputs "behind the scenes" in ArcObjects, so any information that points me in the right direction would be appreciated!


texts: How do I make text Halo in Mapbox Studio work? I tried different width and blur but the text still the same no matter what


texts: Identifying duplicate geometries in different layers with QGIS

texts:
I'm using QGIS 3.16 Hannover under MacOS 10.13 environment. I've got two polygon shapefile layers: 'LIC' and 'ZEPA'. Some polygons (name, code & geometry) are present in both layers and I want to identify them by adding a new field to the layer 'ZEPA' with the value "yes" if that feature is also present in the other layer. I've tried in the field calculator the function overlay_equals(). My try has been
    array_to_string(overlay_equals('LIC', "yes"))

But it returns the new field empty (instead of "yes") for those equal features. What is wrong?


texts: Extracting overlapping polygons using ArcGIS Pro?

texts:
I'm trying to extract polygons on one layer that overlap with polygons on another layer using ArcGIS Pro. I've followed this approach by accessing Related Data.

Indeed, I can view the overlapping polygons, but is there a way to extract these overlapping polygons onto a separate layer?

Can polygons be extracted using the OBJECTID in each related data row?


texts: creating clusters of points with same attributes

texts:
I have a file of some 4000 points, which belong to around 400 different categories (species of rare plants). I would like to remove duplicates within the categories (which could be repeated observations of the same population, at different dates or slightly different coordinates).

Can I perform some kind of cluster analysis but only within each species? I don't want to divide it up into 400 layers! Ideally the limit of each cluster would be around 1000m, and if each point could have the ID of a cluster added as an attribute that would be perfect. Choosing which points to discard from each cluster is difficult as the most recent year would be the best choice but some have the coordinates recorded at better resolution than others, while others contain valuable information in the free text comment attribute.

I'm relatively new to QGIS and probably not able to use Python/ methods involving code without some serious help!

This is an example of the sort of data I'm discussing (working in UK Ordnance Survey Grid Reference) as you see the first three lines are very close to each other but recorded at different times. I want to amalgamate these/ identify them as a cluster (and eventually delete the older records).


texts: Google directions; API returns no results

texts:
I am trying to figure out why my API called to maps.googleapis is returning no results. I inputted my exact same queries into Google Maps (using their interface though, same numbers) and I get directions.

Here is my example:

https://maps.googleapis.com/maps/api/directions/json?origin=43.4833815,-80.539168&destination=43.47242353,-80.54613339&sensor=true&arrival_time=1417650300000&mode=transit

Here is the result from google maps:

https://www.google.ca/maps/dir/'43.4833815,-80.539168'/'43.47242353,-80.54613339'/@43.4780799,-80.5413407,15z/am=t/data=!4m16!4m15!1m3!2m2!1d-80.539168!2d43.4833815!1m5!1m1!1s0x0:0x0!2m2!1d-80.5461334!2d43.4724235!2m3!6e5!7e2!8j1417632300!3e3

Any ideas why this occurring? I am pretty sure my syntax is correct


texts: How to configure settings to alter the default Panels that open on startup?

texts:
I want the "layers" panel to display automatically when I open any project.  I don't want to have to perform View --> Panels --> Check the Box every time I open QGIS.

How can I do this?


texts: Qt-Designer's spinbox does show or store data

texts:
I am new in Qt Designer. I created a user form to show and insert data to a shapefile. Everything works well, except the spinboxes. The fieldtype in the shapefile is integer with fieldlength 10. The shapefile already contains data. In the userform these data in the spinboxes are not shown, the spinbox shows "NULL". When I look in the attribute table data also are missing, even they have been there before. When I insert data the field length is not 10 as it was before, but only 2. The inserted data are not stored.

What's the problem?


texts:
I have a very large shapefile with 130 columns. I'd like to extract this into 130 layers in mapbox studio.

To do this I have split into three files and uploaded them as tilesets. But as a tileset I can't seem to access the columns.

I see there is a Datasets section as well, dont know if I should be uploading as a dataset instead?

Eseentially I would like to create 130 layers which will show building level data. I am going to use Mapbox GL to show/hide layers to show different attributes.


texts: ArcPy from points to lines keeping points attributes?

texts:
I have a script in ArcPy that read a CSV file like this:

FID,LonDD,LatDD,Length,Depth,Type,Label,Date/Time
1,-1.090261,59.887594,,5,A,Text,041100ZAUG19
1,-1.090261,59.854044,,5,A,Text,041100ZAUG19
1,-1.104367,59.831088,,5,A,Text,041100ZAUG19
1,-1.129053,59.785178,,5,A,Text,041100ZAUG19
2,-0.164150843,60.63374438,,1,B,Text,041100ZAUG20
2,-0.380170968,60.52978211,,1,B,Text,041100ZAUG20
2,-0.436680462,60.5024162,,1,B,Text,041100ZAUG20
2,-0.5844909,60.43050015,,1,B,Text,041100ZAUG20
3,0.311851222,60.34566342,,2,C,Text,041100ZAUG21
3,0.031944655,60.31490018,,2,C,Text,041100ZAUG21
3,-0.283586398,60.28072305,,2,C,Text,041100ZAUG21
3,-0.468049486,60.26297731,,2,C,Text,041100ZAUG21

and it creates lines following the FID column (each value of FID correspond to a different line, i.e. FID=1 is the first line and each row is a vertex of the line), but when I open in ArcMap I don't have the attributes.
I lose them when I use arcpy.PointsToLine_management.

Is there a way to keep them? Maybe for each point of the line?

This is an extract of the script I use:

import archook

# Create AML python environment
archook.get_arcpy()
import arcpy

spRef = r"Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj"

arcpy.CreateFileGDB_management(r'D:\out', 'test.gdb')

# Make the XY event and save to layer file
arcpy.MakeXYEventLayer_management(r'D:\in\data.csv', 'LonDD', 'LatDD', r'D:\in\data', spRef, None)

# Save to layer file
arcpy.SaveToLayerFile_management(r'D:\in\data', r'D:\out\line')

# Convert layer file to gdb. file
arcpy.env.workspace = r'D:\out'
arcpy.CopyFeatures_management(r'D:\out\line.lyr', r'D:\out\test.gdb\line_points')

# Determine name of new line feature class - takes name of input file
arcpy.PointsToLine_management(r'D:\out\test.gdb\line_points', r'D:\out\test.gdb\line', 'FID')

arcpy.Delete_management(r'D:\out\test.gdb\line_points')

os.remove(r'D:\out\line.lyr')


texts: ogr2ogr cannot append of truncate PostGIS table

texts:
I am reading in GML data to a specified PostGIS schema.
This is part of an automated job so and the db table might include views so I don't want to drop the table each time the script is run.

EDIT: It append and truncate works if I dont specify the schema. But in my workflow I have to use a specified schema.

Trying to do this:

Typical use case: “ogr2ogr -append PG:dbname=foo abc.shp –config OGR_TRUNCATE YES

ogr2ogr -append -f "PostgreSQL" PG:"host=localhost user=postgres dbname=geospatial password=postgres" \
    -lco SCHEMA=test \
    -lco GEOMETRY_NAME=geom \
    "C:\Users\my_GML.gml" \
    –config OGR_TRUNCATE YES

This result in ERROR 1: Couldn't fetch requested layer 'â€“-config'! ????

EDIT: Double dash is the correct way –-config OGR_TRUNCATE YES

In addition if I run a standard append version

ogr2ogr -append -f "PostgreSQL" PG:"host=localhost user=postgres dbname=geospatial password=postgres" \
    -lco SCHEMA=test \
    -lco GEOMETRY_NAME=geom \
    "C:\Users\myGML.gml"

ERROR 1: Layer test.mylayer already exists, CreateLayer failed.
Use the layer creation option OVERWRITE=YES to replace it.

What am I doing wrong?


texts: Performing object based image analysis of Landsat scene in eCognition Developer?

texts:
I am completely brand new to eCognition Developer and want to test out Object Based Image Analysis (OBIA) on a Landsat scene of one portion of the Chesapeake Bay. I have been unable to find any online resources for beginners to learn the program or how to even begin doing this. Before going any further, is running OBIA even feasible on a medium spatial resolution image like Landsat? If not, I will continue to refine my pixel based classifications.

Note that I only have the trial version but I do have access to all available tools in the program.


texts: Geoexplorer cannot Save and Publish maps

texts:
I installed Geoexplorer 3.0.2 according to instructions.
I also setup GEOEXPLORER_DATA in web.xml file and it seems that it works fine because http://localhost:8080/geoexplorer/maps/? returns:

{"maps":[]}

But when I'm trying to Save my map I get back the half link:

http://localhost:8080/geoexplorer/composer/

Notice there is no #maps/MapID at the end of the above link,as it was supposed to be.

Also when I'm trying to Publish my map I get the following link:

<iframe style="border: none;" height="400" width="600" src="http://localhost:8080/geoexplorer/viewer/#maps/undefined"></iframe>

Notice that instead of MapID at the end of the above link, it gives me the undefined

The directory GEOEXPLORER_DATA and geoexplorer.db have full permissions to tomcat7 user/group.
Is it an insertion problem in my db file? Is it a bug or something else?

If I run GeoExplorer in development mode (without setting any GEOEXPLORER_DATA directory), I have no problem. I can Save and Publish my maps. I only have that problem when I run it under deployment in a production environment (under tomcat7).

Maybe is a GeoServer Proxy problem because during debugging I get:

POST http://myexternal_IP:8080/geoexplorer/maps/ 401 (Unauthorized) GeoExplorer.js:396
    b.send GeoExplorer.js:396
    (anonymous function)
Trouble saving: {"error":"Not authorized"}

My configuration is the following:

I set up Proxy Base URL of Geoserver to http://myexternal_IP/geoserver
I also setup Apache ProxyPass:

ProxyRequests Off
ProxyPreserveHost On
<Proxy *>
Order deny,allow
Allow from all
</Proxy>
ProxyPass /geoserver http://localhost:8080/geoserver
ProxyPassReverse /geoserver http://localhost:8080/geoserver

Finally, I build the app with:

ant -Dapp.proxy.geoserver=http://myexternal_IP/geoserver debug

But the problem still remains...


texts: Why does Spatialite layer prompt for CRS if SRID was not explicitly defined in Layer Options?

texts: In QGIS, if I Save As... a SpatiaLite file and don't enter an SRID in Layer Options, I will be prompted for a CRS upon opening the file. Yet GDAL documentation says that this option is only for forcing a different CRS/SRID, and the driver should normally be able to pick up the CRS and set it automatically. Is this expected behavior for QGIS, or is this a bug that should be reported?