texts:
Background:
 Two years ago I used QGIS to digitize thousands of Alaskan stream segments using  Bing imagery as the basemap, via the Quick Map Services plugin.  The detailed Bing imagery allowed me to trace the streams at a scale of 1:5,000.  The stream CRS was identical to Bing (WGS 84 / Pseudo Mercator, EPSG: 3857).  I was pleased with the results, as was my client.

The problem:  Since then, Bing updated their imagery for a large portion of my project area.  The new imagery is consistently shifted 5 - 6 meters to the east of the previous imagery.

As a result, my digitized streams (still in EPSG: 3857) no longer align with the new imagery, as seen below:

Here's another scene that depicts the mess that I now face.  The client expects the streams to align with the imagery, and I'm going to have to eat the $$$ required to fix this:

Solution:  I would be happy to somehow "roll back" the imagery to the earlier date so that it once again matches my stream lines.   However, previous SE posts (Accessing older Bing imagery in ArcGIS Desktop? and How to avoid the ESRI basemap (Bing imagery) to update?) state that this is not possible.  I now plan to  select the streams within the offset area and apply rubbersheeting and/or manually move them. This is going to be a lot of work, and I'm not happy.

My question:

Before I start this big editing process, how do I know which imagery is spatially correct?

When I initially digitized the streams, I trusted the Bing imagery's location.  But now I don't know what to believe; is the new imagery incorrectly offset, or is the older imagery at fault?

Perhaps both are incorrect!

Most importantly, what happens if (when) the next imagery update does not align; am I doomed to continually "chase" the imagery by applying edits to the stream locations?


texts: How to identify polygons that have shared boundary?

texts:
I have hundreds of polygons that represent plants geographical ranges around the globe.  For each species I have both a shapefile and collection of findings as points.
I would like to analytically find each pair of species that common range boundary at least at one point.  Or, even better, which species has a boundary in certain proximity from another species.
Any idea how to do this?
 I use QGIS and prefer Perl, but solution in any other language would be fine also.


texts: Spatial regression: general nesting spatial model using splm package in R

texts:
I have a shapefile of panel data (multiple years of observations for each polygon). I am estimating fixed effects spatial regression models with this data. I used one year of this data to generate a spatial weights matrix based on the queen's contiguity neighbors:
spMini = dtaSub[ which(dtaSub$Year=='2001'),]
t2 = as(spMini, 'Spatial')
dta = as.data.frame(dtaSub)
neighbors = poly2nb(t2, queen = TRUE)
matrix_list = nb2listw(neighbors)

#Make shapefile into dataframe for regressions
dta = as.data.frame(dtaSub)

I have used the splm package's maximum likelihood-based functions to specify a spatial lag (assumes rho nonzero), spatial error (assumes lambda nonzero), and spatial autocorrelation (assumes rho and lambda nonzero) model. These are as follows

#Design a model
fm = nmPxlsD ~ cmps_13 + cmps_35 + cmps_57 + cmps_79 + cmp_911 +
  cm_1113 + cm_1315 + cm_1517 + cm_1720

#Spatial lag model - estimates lambda
spatialError = spml(fm,
                data=dtaSub,
                listw=matrix_list,
                model="within",
                effect = "twoway",
               spatial.error = "b",
               )

#Spatial error model (SEM) - estimates rho
spatialLag = spml(fm,
                  data=dtaSub,
                  listw=matrix_list,
                  model="within",
                  effect = "twoway",
                  lag = T)

#Spatial autocorrelation model - estimates rho and lambda
spatialAuto = spml(fm,
                  data=dtaSub,
                  listw=matrix_list,
                  index = c("ID", "Year"),
                  model="within",
                  effect = "twoway",
                  lag = T,
                  spatial.error = "b") #"b" = "kkp" for a FE model

GNM = spml(fm,
                   data=dtaSub,
                   listw=matrix_list,
                   listw2=matrix_list,
                   index = c("ID", "Year"),
                   model="within",
                   effect = "twoway",
                   lag = F)

I cannot find any documentation on how to use splm to estimate a General Nesting Spatial Model (which assumes rho, lambda, and gamma are nonzero) using maximum likelihood. Is the function different from "spml"?


texts: Python Scripting in QGIS to add features and select them?

texts:
I'm trying to migrate a python plugin I wrote from QGIS 2.18 to QGIS 3.4. In the plugin, I'm creating new vector features and eventually adding them to an existing vector layer. I'm trying to have the added features selected automatically in QGIS after running the code described below.
For QGIS 2.18, the QgsVectorLayer class had an addFeatures method that contained a very convenient makeSelected flag. When the flag was set to true, the added features were selected.

In QGIS 3.4 the addFeatures method still exists, but the flag is gone.
I also tried to select the features using some of the class's selection methods but with no success. It seems that after adding the features to the layer, their id's change and I don't have access to these specific features and can not select them using their old ids. The only way that somehow works is assigning some sentinel value to one of the fields in all the new features and then selecting them using a query. I'd like to avoid this way because I'm afraid that at some point it can mess up a user's data.

here's an example of how I created features and added them to the current selected layer in Qgis. it can be run in the Qgis Python terminal:

layer = iface.activeLayer()
layer.startEditing()
feat = QgsFeature(layer.fields())
point_list = [QgsPointXY(0, 0), QgsPointXY(0,1), QgsPointXY(1,1), QgsPointXY(1,0), QgsPointXY(0,0)]
feat.setGeometry(QgsGeometry.fromPolygonXY(point_list))
layer.addFeature(feat)

# I'd like to add feat to the layers selection at this point

self.layer.endEditCommand()

I'm using QGIS 3.4.1 (migrating from QGIS 2.18.26), and the pyqgis API of the above versions.

I'm aware that Qgis3 uses python 3 instead of python 2.


texts:
Morgan Herlocker answered in another question, that he uses turf.js on his desktop with node.js.

I'm really new to GIS and won't become pro for all hacks. My need is:

having a geojson with a polygon - I got it
add a buffer of 5km to that geojson-plygon - well, my problem
having again a geojson with only the buffer - the result from 2
add the geojson-buffer to a OpenStreetMap - I got it

I only use buffer and all  operations are the same every time. That is why I do not have to step in deep.

This is what happens after the script-way:

It remains a problem of epsg:3857 or epsg:4326

The result is imported to OpenStreetMap.


texts: how to remove some fields with the WMSGetFeatureInfo function?

texts:
I am displaying the attribute data in the popup on clicking on the map using openlayers,geoext.
But i want to hide some fields(fid,.....) from the table in the popup.
How to do this.
Below is my code

var control = new OpenLayers.Control.WMSGetFeatureInfo({
               url: "http://localhost:8080/geoserver/cite/wms",
                "LAYERS": 'xxxxxx',
                  srs:'EPSG:3857'
           });

        control.events.on({
            getfeatureinfo: function(event) {
                //close existing popup
                if (popup) {
                    popup.destroy();
                }
                popup = new GeoExt.Popup({
                    title: "Feature Info",
                    map: mapPanel.map,
                    lonlat: mapPanel.map.getLonLatFromViewPortPx(event.xy),
                    width: 250,
                    autoScroll: true,
                    collapsible: true,
                    bodyStyle: {padding: 5},
                    location:event.xy,
                    html: event.text
                    });

                popup.show();
                }
        });

 mapPanel.map.addControl(control);
    control.activate();


texts: How can I limit range of the thumbs in esri/dijit/TimeSlider

texts:
From the example, range between those two  thumbs is unlimited.

So, I need some thing like this

timeSlider.setThumbRangeLimit(1);  //one year

or something that limit the range only 1 year.

I have tried to use the single thumb but it didn't work(Look at the Time: field). I think it needs the range for time.

So, do someone has any solution to achieve this?

Thanks in advance :D


texts: Get satellite lon/lat from NMEA sentence


texts: How to project WGS84 to Kertau RSO in ArcGIS Android

texts:
How do I make a conversion from WGS84 (longitude, latitude) to Kertau RSO using ArcGIS Android?

I looking into this for my project. Am I looking at the correct direction?
I have following points in WGS84
103.78627623849339, 1.4367857974669223
and my output points in Kertau RSO
373636336, 5170376

GeometryEngine.project(point p, spatialReference sr1, spatialReference sr2)

Please advise.


texts: Efficiently read large tif raster to a numpy array with GDAL?

texts:
I have been using the Python GDAL API to read tif raster files as NumPy arrays. Previously, I simply read the raster into an array directly with GDAL:

ds = gdal.Open('example.tif')
fullarray = np.array(ds.ReadAsArray())

However, with larger rasters, I receive a MemoryError. As a workaround, I have been looping over the raster and reading windows into a NumPy array:

#arbitrarily choosing rows, cols value here
rows, cols = 20000, 20000
arr = np.zeros((4, rows, cols))
band = ds.GetRasterBand(1)
xsize, ysize = band.XSize, band.YSize
x_edge, y_edge = int(xsize - cols + 1), int(ysize - rows + 1)
x_extra, y_extra = int(x_edge%cols), int(y_edge%rows)

for i in tqdm(range(0, x_edge, cols)):
    for j in range(0, y_edge, rows):

        #read dataset into pre-allocated array
        ds.ReadAsArray(i, j, cols, rows, arr)

I experienced some significant speed-ups (at least 2x faster) after using a pre-allocated array. However, this is still quite a bottleneck in my code. Is there a better way to do this?


texts: ArcGIS python add in to call Python script for importing mdb data

texts:
I'm trying to change from a.net toolbar to a Python toolbar in ArcGIS 10.2 and am having some difficulty.  I have some stand alone python script which I know works and does what I need to do but I can't work out how to get the add in button to call this script or how to incorporate the script in with the button add in script.

import arcpy
import pythonaddins

class ImportCSVFile(object):
    """Implementation for Import_CSV.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        script1 = "C:\Users\xxx\Documents\AddIns\pythonAddin\Install\xxx"
        execfile(script1)

I've also tried using os.startfile(xxxxx) instead of script1 = "xxxx" and execfile(script1).

I don't get an error but the button doesn't do anything.  Does anyone have any ideas of how to get this to work?


texts: Openlayers remove previously drawn polygon on starting to draw next one

texts:
Objective:

At any point of time there SHOULD only be zero/one polygon(s) on a vector layer with a polygon control. I need the previous polygon to be removed as soon as user starts drawing the first vertex of the next polygon.

Code:

var map = new OpenLayers.Map(mapDiv, options);
var tempLayer = new OpenLayers.Layer.Vector("Temp Layer");
map.addLayer(tempLayer);

var polyControl = new OpenLayers.Control.DrawFeature(tempLayer,OpenLayers.Handler.Polygon);

polyControl.events.register('featureadded', ' ', removeExistingPolygon);

function removeExistingPolygon(feature) //Assuming event listener gets a feature object
//Written mainly with featureadded in mind
{
  if(feature.layer.features.length > 1){
    feature.layer.removeFeatures(feature.layer.features[0]);
  }
}

tempLayer.events.register(<event_name>, ' ', <prev_polygon_removal_function>);

Question:

What event_name and function should I use to attain the objective? I have tried 'sketchstarted' on tempLayer but it gets triggered even before user starts to draw a polygon. 'vertexmodified' doesn't seem to work at all.

Objective is partially satisfied when using featureadded with polyControl but in that case the first polygon remains visible till the second one is completed.

Any pointers on how to proceed are most appreciated. Thanks.


texts: Setting multiple WMS layers over basemap of OpenLayers 3?

texts:
I'm developing a web-based application. In this project I'm trying to create multiple WMS image layer over the map with legends.

As you can see below its showing the WMS layer over its basemap .  This is what I am trying to get.

Below I have tried to displayed one sample from GeoServer and second one is for land use/land cover but i'm not able to set basemap for this layer as you can see. There is only WMS layered displayed when I checked on landuse/ landcover button.

Creating one layer over the basemap it is possible to show but I need to set this for multiple WMS Layers.

So my main problem is I'm unable to write code for showing WMS layer over the basemap.  This is my code so far.

<script src="sideicon_file/js/main.js"></script>
<script src="jsuntitled/ol.js"></script>
<script>
    //this is my first layer as you can see from screenshot no.2
var landuseLayer = new ol.layer.Image({
    coordinates: [-124.73142200000001, 24.955967,
                    -66.969849, 49.371735],
    source: new ol.source.ImageWMS({
      url: 'http://localhost:8080/geoserver/topp/wms',
      params: {'LAYERS': 'topp:states'},
     maxZoom: 19
    }),
    visible: false,
  });

  var  soilLayer = new ol.layer.Image({
  coordinates: [79.14511833527447, 20.987418098133496, 79.2018842619151, 21.050233196545],
        source: new ol.source.ImageWMS({
          url: 'http://localhost:8080/geoserver/sagy/wms',
          params: {'LAYERS': 'sagy:pachgaon_LULC'},
          maxZoom: 19
        }),
        visible: false
      });

var map = new ol.Map({
  layers: [landuseLayer, soilLayer],
  renderer: 'canvas',
  target: document.getElementById('map'),
  controls: ol.control.defaults({
                attributionOptions: ({
                    collapsible: false
                })
            }).extend([
                new ol.control.ZoomSlider(),
                new ol.control.ZoomToExtent({
                    extent: [
                        813079.7791264898, 5929220.284081122,
                        848966.9639063801, 5936863.986909639

                    ]
                }),

                new ol.control.Rotate(),
                new ol.control.OverviewMap(),
                new ol.control.ScaleLine(),
                new ol.control.FullScreen(),
                new ol.control.MousePosition({
                    coordinateFormat: ol.coordinate.createStringXY(4),
                    projection: 'EPSG:4326'
                })
            ]),

  interactions: ol.interaction.defaults().extend([
                new ol.interaction.Select({
                    condition: ol.events.condition.mouseMove
                })
            ]),

            view: new ol.View({
                center: new ol.geom.Point([-10997148, 4569099])
                    .transform('EPSG:4326', 'EPSG:4326').getCoordinates(),
                zoom: 4
            }),

            target: 'map'
        });

$('input[type=checkbox]').on('change', function() {
  var layer = {

    landuse: landuseLayer,
    soil: soilLayer,

  }[$(this).attr('id')];
  layer.setVisible(!layer.getVisible());
});


texts: How to build a Query Layer in ArcMap 10.1 using columns of two different databases?

texts:
As you may know, to create a Query Layer in ArcMap 10.1, it is required to specify a connection to a concrete database and then, define the access to the tables and columns in the SQL Query associated.

It turns out that in my query I need to access columns in the Master system database, such as master.dbo.spt_values.number, but once I am connected in ArcMap 10.1 to the database where to perform the query, can I access system databases and use its fields in my query?

Thanks in advance.


texts: Disappearing Field Boundaries in ArcMap?

texts:
I am creating some polygon features that map farm field boundaries. After creating about 3-4 of them, I noticed that the first one I created disappears from the Map View. If I look at the attributes table, it is still there, and its label shows up. I can even select it and the highlighted outline displays. All the other features in the same layer display. I am at a loss.

Any ideas as to why the feature will not display?


texts: Solution to "Error 000918: Cannot retrieve feature class extent"?

texts:
I have an ArcGIS model where I merge two multipoint feature classes, Then I would like to use "Point to Raster" on it. Cell extent is specified through the input DEM which all data derives from. It worked many times, but now after changing the storage location of an intermediate result right before the "Point to Raster" tool, I get the error message "Error 000918: Cannot retrieve feature class extent". What does that mean? I didn't change the input data, but the model doesn't work anymore.

System: Win7, ArcGIS for Desktop Advanced 10.2.1, 3D and Spatial Analyst.

edit: When I run the model to this point (Zusammenführen/Merge (2)) as shown below, the problem does not occur anymore. Why is it not possible to run the whole model at once?

See the specified part of the model below:


texts: Select Features returns There was error executing query in ArcMap when using OLE DB connection?

texts:
I have a spatial enabled OLE DB Provider corresponding to OGIS simple feature specification for OLE/COM revision 1.1

I can add spatial data as ArcMap layers successfully using this provider and by making OLE DB connection in "Add Data" menu in ArcMap TOC.All tools in "Tools" toolbar of ArcMap can work with this layers except "Select Features" tool.

When I want to use this tool to select some features of my data,I receive this error:

"There was an error executing the query"

I found that when ArcMap wants to select features of my data source it requests some properties of provider.

it requests 2 propertysets.I can interpret only one of them and the other isn't documented in MS OLE DB Provider toolkit.The first propertyset Id is

{e9be80e6-ffaa-11d0-9a97080009ec734b}

which is unknown for me and I couldn't find it in "oledb.h",and in this propertyset there is only one property id equal to 0.

The second propertyset is DBPROPSET_DATASOURCEINFO and there are 2 property id for
DBPROP_MULTIPLESTORAGEOBJECTS and DBPROP_STRUCTUREDSTORAGE
and my provider return "VARIANT_TRUE" for first and "DBPROPVAL_SS_ISEQUENTIALSTREAM" for second.

Now,May you help me to find what's wrong with my code?what may be the source of this error.and how can I
find the meaning of this unknown property set?ESRI is a member of OGC and says their products support OGC
standards explicity,so this code can not be a secret code specific to ESRI providers,sure there is something wrong with
my provider but I can't find it.


texts: Applying Definition Query for multiple features using ArcPy?

texts:
I'm trying to apply the following code to extract features that contain these letters 'MW' within the field called 'CODE' for multiple layers. But I'm getting a SyntaxError for line 5 (defintion query).

import arcpy
mxd = arcpy.mapping.MapDocument("Current")
layers = arcpy.mapping.ListLayers(mxd)
for layer in layers:
    if layer.supports(definitionQuery)
         layer.definitionQuery = ""CODE" LIKE "('%MW%')"


texts: Why does GPS positioning require four satellites?

texts:
I have a question on the GPS positioning algorithm. In all books I've read for 3D positioning we need four satellites, and I don't understand why.

We need to calculate three variables: x, y, z. We know when satellite send the signal to earth and when we receive it we can measure the time the signal travel to earth by checking the shift in PRN generator. For what purpose do we need four satellite?


texts: DBeaver spatial output window: Is it possible to sync/filter spatial query to map bounds?

texts:
While upgrading PostGRES/PostGIS from 9.x to 11.x, I have discovered DBeaver and so far I am pleased.

Apparently it has some spatial integration, including the ability to visualize spatial data results using an integrated Leaflet map. Pretty cool. But I'm missing a slick way to support the query with a map bounds filter, or at least something loose like a right-click > "Copy Map Bounds" option on the map object itself, etc. After some experimentation and Googling I'm either not seeing it or I'm not using the right search terms.

Does anyone know if there is a built-in configuration to sync DBeaver's spatial results viewer's map bounds to the active query?

[Edit: 5.8.2019] Response to question in the comments. DPSSpatial asked..

  ..do you mean that if you select * from a table, the results viewer will only show X amount of features, but if you zoom in, it will then show what is in the map view?

Note how, in this second graphic, the software imposes an arbitrary, customizable limit of 1000 records on any potential query--this is the circle in the lower-right. (I think it was originally 200 then I changed it.) That makes sense for either the Grid or Text view states, which are selectable just to the left of the output frame. But for the Spatial view state, to me, it makes more sense to sync the arbitrary limit to the integrated map's bounds, or to at least have that option, which is what led me to ask about this. However based on what I'm seeing in the GitHub link DPSSpatial provided, I'm thinking this feature is so new it just hasn't evolved very far and likely this isn't possible. (Yet?) :)


texts: Setting up ArcGIS Print service for ArcGIS Online account?

texts:
Can anyone tell me is it possible to setup an ArcGIS Print service with ArcGIS Online account? I could not find where to create one.

What I would like is to implement this feature:
https://developers.arcgis.com/javascript/jssamples/widget_print.html

It requires a print service. Something looks like:
http://sampleserver6.arcgisonline.com/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task


texts: QGIS CRS does not have correct coordinates

texts:
I'm using QGIS to extract centroid coordinates in a multi-shape object from a project I was given. I want the coordinates to be regular latitude an longitude, but they're coming out huge (see below).

These locations are in New York City, so x and y should be somewhere around 74 longitude and 41 latitude. I tried to use this:

How to get correct format of coordinates for export to geojson from QGIS?

where I set the CRS layer for my shapefile and centroids to be EPSG:3628
but it doesn't appear to have worked. I have the same results using the CRS found originally in the project folder. I think these are not the problem as adding an OpenStreetMap layer using either CRS displays New York City.


texts: Create new polygon files from existing vector layer in R

texts:
I have an existing vector shapefile with all the counties in the United States (called counties.shp) I also have a dataset that contains many species-county presence data points (species.csv) which looks something like this:
 GEOID  county    state         Species
1: 21007 Ballard Kentucky     plant1
2: 21007 Ballard Kentucky     plant2
3: 21007 Ballard Kentucky     plant3
4: 21007 Ballard Kentucky     plant4
5: 21017 Bourbon Kentucky     plant5

I'd like to create a new polygon for each species that includes all of the counties that they are paired with in the species.csv file. Each row is a unique species-county combination, so both counties and species are repeated throughout that file, yet I'd like to get a single polygon file for each species. There are several thousand species, so I'd like to have a way to do this all at once and compile the resulting polygons in a relatively neat way so I can easily use them for zonal statistics.
I'm imagining this is a somewhat simple task, but I'm completely new to GIS in R.


texts: How to calculate total area of all features of a class in ArcGIS without double counting overlaps

texts: I have a polygon feature class with a field "Area" in the Attribute Table calculating the area of each features. Some features overlap each other. How can I calculate the total area of all the features without double counting the overlaps.


texts:
I use QGIS 2.18.16.

I import 1.000.000 data points, from a text document, as a CSV layer.
Then I add a raster of 5 by 5m over those points.

Project, layer and raster are put in the same CRS (Belgian Lambert 72)
After that I count points in polygon.

This works fine.

How do I count point in polygon for up to 5 million data points or more?
When I import more then 1 million I can't count points in polygon no more.


texts:
I'm trying to create a tool which enables the user to draw a linestring with only 2 points. It's surprising how long I've searched and can't find an example or someone asking the question.

I know that there is a potential maxPoints property but I don't know how to apply it. My guess is that it involves defining a geometry function on the linestring draw interaction but I'm struggling to make it work.

This is what I've got so far:

var vecSource = session.activeRoute.Layer.getSource();
// setup draw interaction
session.draw = new ol.interaction.Draw({
    source: vecSource,
    type: geomtype,
    geometryFunction: geometryFunction,
    maxPoints: maxPoints
});

    if (geomtype === 'Line') {
    maxPoints = 2;
    geomtype = 'LineString';
    geometryFunction = function (coords, geom) {
        geom = new ol.geom.LineString();
        coords = geom.getCoordinates();
        if (coords.length == 2) {
            finishDrawing();
        }
    }
}

I realise this probably doesn't make any sense but I've tried several different methods and this is the best I've come up with. A linestring has been created but length and geom = null when I test this.

Can anyone please offer any insight?


texts: Live KML Feeds - From the Ground Up, How Would You Build It and Why?

texts:
I found the "Live KML from PostGIS to Google Earth" post very informative but a bit daunting for someone new to the open source world of products like PostGIS and GeoServer.

If I am building the product from scratch, is the PostgreSQL / PostGIS - GeoServer combination the accepted best path to take?

What are some alternatives or strengths / weaknesses?


texts: Generating random raster using weights with negative binomial distribution using R?

texts:
I'm trying to generate "species counts" in raster form -- essentially, the number of individuals in a species in a raster cell. I'd like to use the negative binomial distribution to do so. Normally I would just use something like the Create Random Raster tool in ArcGIS Desktop; however, there is an extra wrinkle. I'd also like the species to be weighted based on a second raster, where higher values are more likely to fall in certain areas.

So if I were hypothetically modeling a species like a lion, I'd create a raster with possible values taken from a negative binomial distribution where higher values were more likely on savannahs, less likely in urban areas, and never found in the middle of lakes or oceans (with the land classes taken from an NDVI raster, or something similar).

Is the R raster package capable of doing this?


texts:
I'm trying to use the "Fixed Distance Buffer" from the Processing Toolbox to create a 500m buffer around all the woodlands in my layer, but the tool only puts a buffer around the smaller woodlands.

How can I make it put a buffer around every attribute?


texts:
I have following questions:

  Access Landsat image(s) web services in QGIS
  And can I do analysis based on these WMS?

Is it possible?

if so, when I tried to access from Layers>Add Rasters> the formats shows all but not WMS image service?

Is there any other method to do so?

do I need to install any extension to for WMS images?

or through Python programs?


texts: Removing overlapping areas between two rasters?

texts:
I need some help with mapcalc.

I have two rasters:

A contains cells with 1 or 2 or null;
B has values>=1 or null.

I would like to delete areas where A overlaps with B, then get a resulting map where there is 1 where B has a value>=1 and null where both A and B have >= 1.

It seems easy but I'm trying many different combinations unsuccessfully.


texts: How to place labels in the top left corner of every polygon?

texts: I'm trying to label the polygon boundaries instead of an anchor point. Not sure if this is possible. For this example I wish to label these numbers in the top left corner of every polygon.


texts: Populating layer field with dict.values() in PyQGIS 3?

texts:
I've coded the script below in PyQGIS 3 to get the azimuth of starting or ending points for a line layer that touches a border. It works fine, but I want to populate the azimuth field layer with the value from azDict and doesn't work.

from shapely.geometry import LineString, Point, Polygon, LinearRing

lLayer = QgsProject.instance().mapLayersByName('roads')[0]
pLayer = QgsProject.instance().mapLayersByName('aux_point')[0]
polLayer = QgsProject.instance().mapLayersByName('aux_boundary')[0]
feat = QgsFeature()
pr = pLayer.dataProvider()
vDict = {}
azDict = {}

# iterates in vertex lines to get dict = {Point1:Point2}
for road in lLayer.getFeatures():
    if road.geometry() is not None:
        line = road.geometry().asMultiPolyline()
        for vrtx in line:
            stp = QgsPointXY(line[0][0]) #  start point
            scp = QgsPointXY(line[0][1]) #  second point
            vDict[stp] = vDict.setdefault(stp, scp)
            edp = QgsPointXY(line[-1][-1]) #  end point
            lbo = QgsPointXY(line[-1][-2]) #  last but one point
            vDict[edp] = vDict.setdefault(edp, lbo)

#  iterates in vDict to calculates the azimuth
for k,v in vDict.items():
    azim = k.azimuth(v)
    azDict[k] = azDict.setdefault(k,azim) #  gives a dict = {Point:Azimuth}

#  begin the spatial filter
for polygon in polLayer.getFeatures(): #  transforms my boundary in a shapely LinearRing
    if polygon.geometry() is not None:
        poly = polygon.geometry().asMultiPolygon()
pol = Polygon(poly[0][0]) #  creates a shapely polygon
linearring = LinearRing(list(pol.exterior.coords))

#  iterates in azDict to filter spatially
for k,v in azDict.items():
    if linearring.distance(Point(k)) < 1:
        orientation = pr.fieldNameIndex("azimuth")
        feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(Point(k).x,Point(k).y)))
        pr.addFeatures([feat])
        #  try to write the azimuth in the field but doesn't works
        pLayer.startEditing()
        pLayer.changeAttributeValue(feat.id(), orientation, v)
    QgsProject.instance().addMapLayer(pLayer)
    pLayer.triggerRepaint()
    pLayer.commitChanges()

Am I losing something?


texts: Choosing strategic spatial data model for Electrical Utility sector?

texts:
Is there an ideal industry standard spatial data model for an electrical utility company (and which utilities use it) that has a mix of linear & point assets or is it better to develop a custom model?

What are the pros and cons for adopting an industry standard model in this area?


texts: Why isn't the QGIS 3.4 tracing tool not tracing my line?

texts:
I'm trying to trace a new line onto another existing line layer (highlighted yellow), but the red rubber band will not stick or follow the existing line completely. Instead the red rubber band wanders around to points on the existing line layer or other layers nowhere near where I'm trying to trace. The map project was created in QGIS 2.18. Snapping on/all layers/vertex/tracing enabled, no offset. Snap tolerance 10px.
I tried this on a short simple line in another project and it worked, but not on the project layers I'm working with. Picture shows the chaos. AutoTrace plugin worked fine in 2.18.
Here's a pic to see how the rubber band is behaving.


texts: QGIS v. ArcMap Network Analysis Isochrones

texts:
I am looking at developing catchment areas based on a street-lines feature class against various bus stops to find an ideal catchment area which can then be converted and run in a simulation. The impedance is distance over time such that a vehicle will only make it so far, by x-time. The images below show a 5,10,15-minute area. The visual aesthetic of the ArcMap representation is a bit shoddy but I just wanted to have something for comparison.

In ArcMap I am able to control the speed attribute for the network impedance using attribute table values in the Network Analysis extension. I need to somehow control QGIS interaction with the road network to 25 mph if that's possible.

QGIS Plugin ORS Tools (previously OSM Tools) only allows for an input classification such as 'driving-car' in the Isochrones generation filed under the Batch Jobs tab. I need to limit the speeds of the road network to 25 mph and I have not found a good way of doing this.

Comparisons shown below are using the same point in the City of San Jose, California. The QGIS image shows the point in red. I forgot to make the point visible in ArcMap. The ArcMap image is different because it is able to control for the street-speeds.

I have done several forms of research on this topic and ORS was the one tool/plugin that seemed to accomplish what I was looking for.

How can I set the speed to 25 mph using QGIS?

ArcMap

QGIS with ORS


texts: Load layer in QGIS legend from Processing script [QGIS 3]

texts:
I have a weird problem. I'm writing a Processing standalone plugin in QGIS 3 (current master, 3.1).

In the main algorithm I have to create some intermediate layers (mostly raster) with GRASS and SAGA provider, no problem so far. For example here the snippet of r.watershed:

viewshed = processing.run('grass7:r.viewshed', {
    'input': raster,
    'coordinates': coordStr,
    'observer_elevation': 0,
    'target_elevation': 0,
    'max_distance': Visibility_Distance,
    'refraction_coeff': 0.14286,
    'memory': 500,
    '-c':False,
    '-r':False,
    '-b':False,
    '-e':False,
    'output':os.path.join(Results, 'my_output.tif'),
    'GRASS_REGION_PARAMETER': tr_extent,
    'GRASS_REGION_CELLSIZE_PARAMETER':0,
    'GRASS_RASTER_FORMAT_OPT':'',
    'GRASS_RASTER_FORMAT_META':''

})

but then I cannot load the resulting layer in the legend, even if it is correctly created and I can load it by hand.

Here what I'm trying to do:

view = QgsRasterLayer(viewshed['output'], 'my_output','gdal')
QgsProject.instance().addMapLayer(view)

I don't get any error from the console.

The exact same syntax works in a single script written in the Python editor of QGIS.

Any clue?


texts:
My primary computer has a small screen, and QGIS keeps coming up with 4 rows of toolbars, at least 3 of which have large empty space at the end.
Currently the last row has only one icon.   Is there some way I can move this icon into some of the empty space in previous rows?


texts: Dragging all vertices at a given point programmatically in PyQGIS

texts:
I'm writing a Plugin where you build networks with point and line geometries. One Tool should give the User the possibility to make changes to the network by dragging and dropping the points to another location. Since the points are sometimes connected to lines I want the vertex of the connected lines to the new location as well. I managed to write a tool that does exactly that, but if several line features are connected to the same point it fails and only drags one line. I tried looping through all the line features and thought it should work, but it still only moves one line and I don't know why.

Here a few pictures to illustrate my problem:

This is how the network could be made up

changing one point with one line works as expected

changing the middle point fails

I wrote a maptool-class passing the QgsMapTool to it.

In the canvasPressEvent I call on two methods one looks for features and the other for vertices. if there are results i enable dragging and start editing. Here are the code parts where I think the problem might be in

def findVertexAt (self, pos):

    self.lineFeatures=[]
    self.vertices = []

    point = QgsPoint(self.canvas.getCoordinateTransform().toMapCoordinates(pos.x(), pos.y()))
    tolerance = self.calcTolerance(pos)
    searchRect = QgsRectangle(point.x() - tolerance,
                              point.y() - tolerance,
                              point.x() + tolerance,
                              point.y() + tolerance)
    request = QgsFeatureRequest()
    request.setFilterRect(searchRect)
    request.setFlags(QgsFeatureRequest.ExactIntersect)

    for feature in self.netLayer.getFeatures(request):
        self.lineFeatures.append(feature)
        i = 0
        while i <= len(self.lineFeatures)-1:
            vertexCoord, vertex, prevVertex, nextVertex, distSquared = self.lineFeatures[i].geometry().closestVertex(QgsPointXY(point))
            self.vertices.append(vertex)
            i = i + 1
        return self.vertices

def canvasMoveEvent(self, event):

    point = QgsPoint(self.canvas.getCoordinateTransform().toMapCoordinates(event.pos()))
    pointXY = QgsPointXY(point)

    if self.dragging:

        geometry = QgsGeometry.fromPointXY(pointXY)
        self.layer.changeGeometry(self.feature.id(),geometry)

        if self.lineDragging:
            i = 0
            while i <= len(self.lineFeatures)-1:
                lineGeometry = self.lineFeatures[i].geometry()
                lineGeometry.moveVertex(point, self.vertices[i])
                self.netLayer.changeGeometry(self.lineFeatures[i].id(), lineGeometry)
                i = i + 1

The canvasRelease event then just disables dragging and commits the changes


texts: Calculating Shannon index of habitat diversity

texts:
I would like to calculate Shannon's diversity index for habitat diversity, i.e. the diversity of landscape factors. I have 11 locations (paired fields) for which I already defined (with QGIS) and calculated the proportions of 5 landscape factors (crops, grassland, semi-natural, forest and others) within a certain radius. I would like to compare the diversity of these five landscape factors among the 11 sites, but I have never done this before. I know how to use the R package "vegan" for calculating diversity indices with species, but not with landscape factors.

Could somebody please explain me how to deal with this?


texts:
Can somebody please tell me how to download data from specific locations on the sentinel-1 data hub?

I cannot seem to get the geospatial search function to work...


texts: R: How to standardise extent of rasters for raster stack

texts:
So I'm trying to build a raster stack of 7 variables. In ArcMap I've already standardised resolution and CRS, and I tried to standardise extent and column/row number by extracting by mask with a polygon of the study area. Despite this, extent and column/row number still varies.

Originally I skipped over this issue using quick=TRUE when building the raster stack, but since doing that will create more problems later I decided to try to sort it out now. I've tried using the solution posted here, and ended up with this code:

var_files <- list.files(file.path("D:/OneDrive/Diss/final_layers/"),
                        pattern="tif$", full.names = TRUE)
r.var <- list()
for(i in var_files) {
  ras <- try( raster(i) )
  r.var[i] <- data.frame(name = i, nrow=nrow(ras), ncol=ncol(ras), res=res(ras)[1],
                         proj=proj4string(ras), xmin=extent(ras)[1],
                         xmax=extent(ras)[2], ymin=extent(ras)[3],
                         ymax=extent(ras)[4])
}
do.call("rbind", var_files)
var_stack <- stack(var_files)
plot(var_stack)

And then I get these warnings after the loop:

In addition: Warning messages:
1: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length
2: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length
3: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length
4: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length
5: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length
6: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length
7: In r.var[i] <- data.frame(name = i, nrow = nrow(ras), ncol = ncol(ras),  :
  number of items to replace is not a multiple of replacement length

And this error after do.call:

Error in do.call("rbind", var_files) : second argument must be a list


texts:
I am trying to make a depth map of one of the lakes that I manage. I recorded my data using a lowrance depth finder. I have imported the data into QGIS and now have all the points that were recorded with depth and Lat. and Long. When I try to use the contour plugin I get an exception struck #### error. The error displays a number but when I delete that line from the attribute table the number changes. I keep deleting the number of the error that is returned but it seems like an endless loop.

What am I doing?

Data:

Error:

I can go through and manually delete points and after deleting a certain amount I no longer get the exception struck error. But by the time I do this I have lost many points and the map will no longer be accurate.


texts: How to find out the length and friction costs of my least cost path?

texts:
I modeled different possibilities for transmission lines with the GIS tool "Cost Path". Now I would like to find out more information about my results. First of all: Where can I see the length of my least cost path? Second: I generated my input costs raster based on different weights for different land covers. Now I would like to find out how many km my transmission line goes through each land cover type. How can I do that? (In my weighted overlay, I gave every land cover type another weight).

I am using ArcGIS 10.2.2


texts:
My data has a date field formatted as "MM/DD/YYY hh:mm:ss". I want to create a new field with just the time in the full "hh:mm:ss" format. I've found how to extract only the hour, only the minute, or only the second, but how do I extract the full time stamp?

I've tried using the Field Calculator function DatePart("hh:mm:ss",[DateField]), but to no avail.

Also, should this new field showing only time be a date field, or would it be a string?

Per the possible duplicate question edit; the answer given was to use the "Convert Time Field" tool in ArcMap. However, this tool does not give an Output Time Format option of hh:mm:ss


texts: Using multiple basemaps with masking/clipping?

texts:
Is there a way to limit a basemap to within a boundary in arcmap 10?

I am working on a project within San Joaquin county in California.  I want the areas surrounding the county to be displayed with a light grey basemap with minimal details.  Within the county itself, however, I want to use a more detailed basemap that has the information I am trying to portray.

Is there a way to add a second basemap but have it be limited to the county boundaries?


texts:
I am trying to parse csv file and visualize it as geotiff using gdal_grid but I keep getting this error:

C:\work\testgdal1>gdal_grid -ot Float32 -l test  test.vrt dona.tif
ERROR 1: Failed to open datasource `test.csv'.
Unable to open input datasource "test.vrt".
Failed to open datasource `test.csv'.

Both test.csv and test.vrt are in the same folder.

What can possibly go wrong and how to fix it?

Environment:

Windows 7 x64
osgeo4w as of 2 May 2011

test.vrt

<OGRVRTDataSource>
    <OGRVRTLayer name="test">
    <SrcDataSource>test.csv</SrcDataSource>
    <SrcLayer>test</SrcLayer>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="field_10" y="field_9" z="field_20"/>
    <GeometryType>wkbPoint</GeometryType>
</OGRVRTLayer>

2 lines of test.csv

200010  207  020311    40658.5  406593 52 344927.31 7100203.50  -26.2078720  127.4491855 345060.64 7100369.14   26.4  650.3  628.0 55471.293    20.168 55648.817 55637.523  -146.062
200010  207  020311    40658.6  406594 52 344932.31 7100203.50  -26.2078726  127.4492355 345065.64 7100369.14   27.2  650.3  627.2 55456.719    20.172 55648.814 55637.520  -160.629


texts: Issues with building a VRT. from OS map tiles in QGIS

texts:
I am trying to speed up rendering of OS raster map tiles in QGIS.

The tiles are 1:10k Tiffs - about 130 of them.

When I attempt to build the VRT it says complete, then it shows up an error message:

Warning 6: gdalbuildvrt does not support rotated geo transforms.

The only settings I had on whilst creating the VRT were: resolution - highest
                                                         Add to canvas after finishing

As you can see from the image some of the tiles were successful. I didn't have a problem with creating a VRT for 1:50k or 1:250k map tiles.

Any suggestions or advice on speeding up rendering of OS map tiles (such as pyramids) would be much appreciated.

Thanks!


texts:
I am trying to find an easier method to align the symbols from TASVEG 3.0 with those that I have in my layer in QGIS. I can do it manually and alter each individual symbol that I have categorised but there are over 150 vegetation communities (different symbols).

I am wondering if anyone has an easier way to do this or knows of a download file possibly an SVG with the symbols preloaded.


texts: How can I define the layer CRS for a memory layer in pyqgis?

texts:
How can I define the CRS for a memory layer from the QGIS project settings?

    vpoi = QgsVectorLayer("Point", "vectpoi", "memory")
    feature = QgsFeature()
    feature.setGeometry(QgsGeometry.fromPoint(QgsPoint(10,10)))
    provider = vpoly.dataProvider()
    ...

And how can I define the CRS from another loaded layer?


texts:
I am using the rLiDAR package to read some LAS files, but the readLAS function in rLiDAR is giving the following error:
Error in readLAS("D:/LAS_Files/test_1.las", short = TRUE) :
  no records left to read

However, the function readLAS works just fine with test.las (small file-162 MB), but throws an error with the test_1 (big size, >500 MB)
.The code is as follows:
#Read the LAS files
library(rLiDAR)
las<-readLAS("D:/Summer_2019/NRS_Project/01_data/01_las_2013/test.las",short = TRUE)
las_1<-readLAS("D:/LAS_Files/test_1.las",short = TRUE)


texts: All Network from fixed distance of a point

texts:
I am new to Network Analyst.

I have a network dataset and I want to generate all polylines (streets) calculated from a single origin (point) with a fixed distance (non-euclidean).


texts: I have just downloaded the newest version of QGIS (3.2) and the processing toolbox is nowhere to be found. It says that it is now apart of the core package, but i cannot find it on anywhere on the program.


texts: Solar radiation calculation from GTOP030 DEM image

texts:
Can you help me to calculate solar radiation from GTOPO30 DEM image at the 30 m spatial resolution.

I've tried to do this many times using spatial analysis solar radiation calculation tool in ArcGIS. However, there is always an error message saying unknown problem.

What does it mean? Is there anyone who has a experience with this. Your suggestions will help me a lot.


texts: I have a raster file that show forest cover, very simple, where there is a green pixel there is forest. I need to count how many of this pixels fall under a certain elevation. I have the elevation already calculated and so on. I only need to relate these two features. Like, how much forest are there in different elevations. I am just starting working with QGIS, so I am sorry if the questions is too stupid.


texts: Strange Pattern in LIDAR derived DSM raster

texts:
I produced a DSM raster from LIDAR data for eastern Kentucky. There are strange lines from southwest to northeast. Do you have any idea what those might be?
I generated the DSM raster with both first return and non ground classes but results were similar. I am using ArcMap 10.6. Basically I used first return or non ground classes and created DSM from elevation data by using LAS dataset to raster tool. Then subtracted DEM from it. DEM looks fine. I also created Canopy Density Model (CDM) by using LAS Point statistics as raster (points above surface/ all points). I still have the same issue with CDM. Red areas show zero value (ground) green areas are higher areas with structures or canopy.
Here is the picture


texts: QGIS Qfield visiblity of offline layers on Android 6.0 Samsung tablet Galaxy tab S2

texts:
I am running a brand new out of the box Samsung Galaxy Tab S2 with Android 6 and trying to get QGIS QField to operate correctly. So far I have created a QGIS project in QGIS desktop, packaged it with QfieldSync and uploaded it to the tablet, but I'm struggling to do anything with it.

It opens and I can see the individual layers in 'Digitise' but only the layer that I copied is visible in the data frame, the layers that I wish to edit and that are in 'offline editing' configuration are completely invisible. Is this normal (I am used to ArcGIS collector - where you can edit data and see it at the same time) ? I'd like to be able to see the existing data in these layers while I edit it. Is this possible and have a made an error or is there some other way to do this?


texts: How to label features in QGIS using SQL statement?

texts:
In QGIS, I can't seem to find a way to define the features that I want to label using a SQL query?  Is this possible?

For example, I downloaded the OSM Natural shapefiles, and a lot of those features have no values in their "name" field.  I would like to display all features, but only label those features where NAME <> "".  The same concept could be applied to roads, where TYPE = "Motorway".

In the image below, I want all polygons visible, but don't want to label the 3 polys on the right.


texts: WFS Query in QGIS to load features within 1km of points stored locally

texts:
I have a points layer saved locally containing sites of interest. I would like to be able to load features from a WFS source that are within 1km of my local layer. I'm looking for something similar to the following in the build query option of the add WFS layer dialog:

SELECT wfs.* FROM wfs, local_points
WHERE ST_DWithin(wfs.geometry, local.geometry, 1000)

However, this throws an error that local_points is unknown. I'm guessing that is because QGIS is own looking for tables on the WFS. Can it look for local tables? Is there any way of achieving what I want? Requesting features that overlap the current extents isn't suitable as the sites are a long distance apart and I'd end up downloading most of the feature set from the WFS.


texts: How to populate and add a new column in CSV file for QGIS Python plugins

texts:
I am trying to add a new column(field) in a CSV file and also populate the CSV file with the result of the multiplication between the 4th column(field) and 0.3048. my code copies everything in the CSV but failed to perform the operation in the line==>: ### append new column with feature attributes.
My plugin code:
filename = self.dlg.lineEdit.text()
with open(filename, ‘w’) as output_file:
selectedLayerIndex = self.dlg.comboBox.currentIndex()
selectedLayer = layers[selectedLayerIndex].layer()
fieldnames = [field.name() for field in selectedLayer.fields()]
# write header
line = ‘,’.join(name for name in fieldnames) + ‘\n’
output_file.write(line)
## wirte feature attributes
for f in selectedLayer.getFeatures():
line = ‘,’.join(str(f[name]) for name in fieldnames) + ‘\n’
output_file.write(line)
### append new column with feature attributes
open(filename, ‘w’) as output_file:
writer = csv.writer(output_file, lineterminator=’\n’)
reader = csv.reader(input_file)

all = []
row = next(reader)
row.append(‘Tile(Elev)’)
all.append(row)

for row in reader:
row.append(row[3]*0.3048)
all.append(row)

writer.writerows(all)


texts: how to add data from a layer to postgis table with different field names

texts:
I am importing waypoints from GPS (in gpx format) using Layer>Add Vector layer and then I want to clean up the data and add some more attributes and finally add them to a master table in PostGIS.

The import works fine but I get a table with a lot of fields I don't need and the fields I do need in some cases have the wrong names nor can I modify the table without doing a save as and reloading it.

So I end up

importing the layer,
saving it as a in a format that I can modify
importing it again
clean up the data and add new attributes (using some of the fields i don't need)

Now I want to update the master table.  The only way I can see of doing this is to save the layer as PostGIS (or do this in step 2) and then use SQL to do the update so I can select which fields I want inserted in which columns.

Is there a more stream lined way of doing this?

PS.  I am guessing the answer is write a script to do it.  To this end I have just purchased The PostGIS cookbook and the Pyqgis book but it will be a few weeks while they travel around the world.


texts:
I need to custom icon for my app . One icon for a group of markers. I add a URL for each icone and I don't know how to add it dynamically to my iconUrl. It to change for each select.

Here is my code :  I tried this but I don't even know if it works!

function displayIconByCAt(iconUrl) {
    console.log(iconUrl);

 var data = iconUrl.toString(data);
    console.log(data);
    createNewIcon(data);
}
function createNewIcon(data){
    var iconMarkers = L.icon({
        iconUrl:'',
        iconSize: [38, 38], // size of the icon

        iconAnchor: [22, 94], // point of the icon which will correspond to marker's location

        popupAnchor: [-3, -76] // point from which the popup should open relative to the iconAncho
    })
    var iconNew = iconUrl;
    var value = data;

    console.log(iconNew);
    console.log(value);

    iconMarkers [iconNew] = value;

    console.log(iconMarkers[iconNew]);

}
console.log(iconMarkers);


texts:
I'm developing an iOS app for displaying my own custom map data (on GeoServer). The SDK I'm using to display my map is route-me (Alpstein fork) as the map viewer for now. route-me supports raster tiles (slippy tiles) and I'm OK with that for the most part, mostly for online mode. But I want to support both online maps and offline routing, plus, caching z17+ is a very bad practice because of the growth in number of tiles and the fact that very little part of these tiles will be viewed by users, and z17+ doesn't necessarily mean new data! we can use the same data and over zoom it.
As part of my research, I have found that OSM protobuf format is a very efficient format for offline usage and routing but I cannot find any good renderer for the format!
My scenario is: for z0-z17 I will load tiles from server with no problem ->from z17-z20 I will try and load a larger set (bbox from z15 of the tile, meaning 16 tiles around) and try to render the tile's bbox from osm.pbf and pass it to route-me as the tile and cache it like a normal tile.

The question now is: How can I render osm files in iOS?


texts: Fire click event on Leaflet layer using latlng pair

texts:
I have a button on a leaflet map that when clicked zooms the map to the users current location. This function is working properly.

<button type="button" id="locator" class="btn btn-primary btn-lg">
..
$("#locator").click(function() {
  map.locate({
    setView: true,
    maxZoom: 16
  });

Using the lat and long of the found location I would like to emulate a click event on the intersecting polygon layer.  To perform this I use the maps 'locationfound' event to fire a click event at the found location (below).  Nothing happens.

map.on('locationfound', function(e) {
  var latlngPoint = new L.LatLng(e.latlng.lng, e.latlng.lat);
  map.fireEvent('click', {
    latlng: latlngPoint
  });
});

The polygon layer has the following click event.
The layer click event works fine when clicked manually, but does not execute via passing latlng pair.

var sections = L.geoJson(randall, {
  style: featureStyle,
  onEachFeature: onEachFeature
});

function onEachFeature(feature, featureLayer) {

  featureLayer.bindPopup(feature.properties.ANUM);

  featureLayer.on('click', function(e) {

    console.log('Layer clicked!', e);
  })
};

Heres my jsfiddle.


texts: Exception for error produced while running arcpy.MakeXYEventLayer_management method

texts:
I am using Python 3.6.8 to make an XY Event layer. I first use the arcpy.MakeXYEventLayer_management(inputCSV, "Longitude", "Latitude", TempFeatureLayer, WGS84, "")

where TempFeatureLayer is defined earlier as "XYTemp", and WGS84 is the variable holding the GCS_WGS84 coordinate system info.

I pass this and get an error the first time that is pretty extensive but the end says: "

  SystemError:  returned a result with an
  error set"

I am able to continue running the code successfully inspite of this seeming error. I would like for this to simply not show up as the code will go to someone with little experience and I dont want this to trigger concern.

Interestingly when I run the code again for the next input, this error does not appear again....

Does anyone know of an exception handling i can use to eliminate this?

I am currently running 10.7 desktop and have been working in the jupityer notebooks version that comes with this desktop version.

I have attached the error and the Python code I wrote. the section happens at line 277

workingDir = r'C:\Users\e17864\Desktop\CPSE_Analysis\Datasets\ByYear\Incidents_'+year#Change
workingGDB = os.path.join(workingDir, "CPSE_FinalDatasets.gdb")

ap.env.overwriteOutput = True
ap.env.workspace = r"C:\Users\e17864\Desktop\CPSE_Analysis\Datasets\CPSE_FinalDatasets.gdb"

WGS84 = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],\
               PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];\
               -400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;\
               0.001;0.001;IsHighPrecision"
XYEventsTemp = "XYEvents" + year #Change Name
XYEvents = "XY_WGS" + year #Change this name
ProjectedLayer = "IncidentData" + year + "_NAD83" #Initial Dataset projected into NAD83_Harn_NM State Plane

XYEventsFL = ap.MakeXYEventLayer_management(CleanedIncidentsPath, "Longitude", "Latitude",XYEventsTemp, WGS84, "")

#Makes XY Event Layer real- prep to join district and parcel data
ap.CopyFeatures_management(XYEventsFL, XYEvents, "", "0", "0", "0")

Error:
SystemError                               Traceback (most recent call last)
C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\pandas\core\ops.py in f(self, other)
   2089     def f(self, other):
   2090
-> 2091         other = _align_method_FRAME(self, other, axis=None)
   2092
   2093         if isinstance(other, ABCDataFrame):

C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\pandas\core\ops.py in _align_method_FRAME(left, right, axis)
   1949         return right
   1950
-> 1951     if isinstance(right, np.ndarray):
   1952
   1953         if right.ndim == 1:

SystemError: <built-in function isinstance> returned a result with an error set


texts: Intersection between two SpatialPolygonsDataFrame using gIntersection function in R

texts:
I want to select a subset of my map (a SpatialPolygonsDataFrame) to create a new vector map = part of the initial map included inside a given polygon (eventually, inside several polygons).

I have a large study area, under the object "parcel", and would like to select the part of the map inside the polygon "hr_kud". This is basically the part of my study area used by an animal, so I want to compute basic statistics on this new map and used it for new analyses (of habitat selection).

I have to repeat this operation for about 100 animals...

Here are my data :

http://dl.free.fr/vX4TA8nsE

parcel = my study area map. and hr_kud = an example of homerange used to reduce the study area

Here is a part of my R code, used under Rstudio (R version 2.15.2 (2012-10-26)) :

  library(rgeos)

  load("Data.RData")

  class(parcel)

    [1] "SpatialPolygonsDataFrame"

    attr(,"package")

    [1] "sp"

  class(hr_kud)

    [1] "SpatialPolygonsDataFrame"

    attr(,"package")

    [1] "sp"

  object.size(hr_kud)

    9536 bytes

  object.size(parcel)

    88090144 bytes

  plot(parcel, axes=TRUE)

  plot(hr_kud, col="red", add=TRUE)

And when I run my final command, it causes Rstudio to crash :

  hr_map <- gIntersection(parcel, hr_kud, byid=TRUE)

It is not the first crash reported, but non of the answers find on the net are efficient for me... It's exactly the same using R and all my packages are updated.

So If you have any idea of why it's crashing, or of an alternative method to create my subset map or intersect two SpatialPolygonsDataFrame in R, I will be really appreciated !

thanks in advance,

Sophie


texts: Using ArcPy to get access to connection details of layer connected via Direct Connect?

texts:
I am working on a python script that loops through each layer in a map document and prints out the database connection details for the particular layer.

I have some layers on the map which require a user/password for the layer to load.

What is the best way to get access to the connection details e.g. Instance, Database platform etc  of a layer that requires a password.

I can get access to the datasource using layer.dataSource without a password but when I try and use arcpy.describe(layer.dataSource) I am prompted with the dialog to enter the username/password.
I would like to get the highlighted information before the dialog appears.

Does anyone know the best way to do this?


texts: Attribute table empty in QGis, although google earth displays info

texts:
I am visualizing a .kmz in google earth and when a point is clicked one can read the attributes, such as lat/long, name, etc.

However, when I open the same file in QGIS the point display, BUT all attributes in the attribute table show as "null", and the original attribute categories are missing.

How can make the attribute table show correctly in QGIS?


texts: Export a large raster to a comma or tab delimited text file in ArcGIS 10.2 for Desktop

texts:
I have a large raster with multiple attributes in 200-meter cells, and I would like to export this to a tab or comma delimited file which contains the X,Y coordinates and all of the attributes.  If it's possible I'd also like to have the file arranged so it goes East to West across the whole raster like a typewriter.

I'm currently in the process of trying raster to point,  assigning X,Y, to the points and exporting the table.


texts:
The South African National Mapping Agency acquires multispectral imagery (R,G,B and NIR)and currently compresses its aerial imagery with JPEG compression with a TIF wrapper. So we give out JPEG (Q3) compressed TIF files. This has a serious limitation in that it does not support 4 band compression. Due to this limitation we give our RGB.tif as well as CIR.tif. This situation has obvious duplication of bands.

We are thinking of moving to 4 band JPEG2000 compression for our workflow as well as giving clients 4 band orthorectified jp2 imagery.

Does JPEG2000 support this 4 band combination of (R,G,B,NIR) ?


texts: Applying MODIS Pixel Reliability and Vegetation Index (VI) Quality mask for MOD13Q1 using ArcGIS Desktop?

texts:
I have some sets of MODIS MOD13Q1 data sets during March and December for year 2000 to 2015. I need to extract the NDVI and EVI for those two months from the entire 15 years.

How do I apply the Pixel reliability and VI Quality to mask out pixel that are unreliable using ArcGIS 10?


texts: Accessing composite data type parameter when validating parameters in Python Toolbox?

texts:
I have a simple tool with two parameters:

a boolean (just to tick or untick it)
a composite data (where I can choose an output location - either a folder or a feature dataset).

I am unable to access this composite data type parameter correctly. I want the script to work like this - if I tick the boolean, I´d like both choices of the composite parameter to be disabled (=grey). If I untick it, I want the second parameter to be enabled.

The code looks like this:

class Tool(object):
   def __init__(self):
      self.label = "Tool"
      self.description = "Prints output location"
      self.canRunInBackground = False

   def getParameterInfo(self):
      param0 = parameter("a", "Tick", "GPBoolean", "Optional")
      param1 = parameter("b", "Output location", ["DEFolder", "DEFeatureDataset"], "Optional")

      params = [param0, param1]
      return params

   def updateParameters(self, parameters):
      if parameters[0].value:
          parameters[1].enabled = 0
      else:
          parameters[1].enabled = 1

      return

   def execute(self, parameters, messages):
      a = parameters[0].value
      b = parameters[1].valueAsText

      if a:
          arcpy.AddMessage("A is ticked.")
      else:
          arcpy.AddMessage("B is chosen instead:" + b)

      del a
      del b

      return

The code in updateParameters is not working as expected, whether I have the boolean ticked or not, the second parameter is still enabled.

I tried parameters[1][0].enabled = 0 but validation says that the object is not subscriptable. Then I tried changing the second (composite data type) parameter to a non-composite data type parameter and it worked. But when it is a composite data type parameter, the parameters[1].enabled = 0 within updateParameters function does not work.

How to make the second parameter disabled when the first parameter is ticked?


texts: Calculations from Feature Server and display values

texts: Using the ESRI JS API I would like to apply some basic math based on the Shape.area of a feature and then return that calculation to the web-page.  I can't see any reference to something even remotely similar to this on the ESRI API Reference page, so I am thinking that this will have to be done entirely using DOJO or DIJIT.  Can somebody please just nudge me in the correct direction.


texts: How to extract depth data of nodes of unstructured mesh in QGIS?

texts:
I am currently working on a project which requires unstructured meshing of a given area. Also, i need to extract the data such as nodes x-y coordinates, cell data(triangular cells in my case) and depth at each node. I am using QGIS for this. I was able to generate mesh through the QGIS-GMSH plugin. I am also able to extract the nodes and cell data but I am unable to extract the depth data.

I have googled a lot to find a solution, but till now, no luck.


texts:
I have MULTIPOLYGON in my PostGIS database to store countries (with islands etc.). But I need to add holes to my borders for lakes. How can I store multipolygon with holes?

I am using ST_SetSRID(ST_Multi(ST_GeomFromText('POLYGON(...)')), 4326) to create multipolygon


texts: Filter Layerlist in Snapping Options Dialog

texts:
I have some large QGIS-Projects where I would like to filter the layerlist that is used for the advanced snapping options.
Even if there are a lot of layers used in the project I only have to adjust the snapping options for some of these layers.

I know that I could set the snapping-options programatically for single layers but I need to adjust the snapping options dynamically during my work (snapping mode, tolerance and so on) so this is not an option.

I think I found the item that is used to display the layerlist:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *

dlg = iface.mainWindow().findChild(QDialog, 'SnappingOption')

st_widget=dlg.findChildren(QStackedWidget)[0]

test= st_widget.children()[2].children()

l_widget= test[1]
l_widget.setVisible(False)
print l_widget

But how can I filter this list in order to show only some of my layers?


texts: Torque time legend incorrect compared to dataset timestamp

texts:
I am creating a map within the CartoDB map creator based on a dataset (spatially located within CT timezone) with a timestamp range of "2015-11-19 11:00:00' AND '2015-11-19 17:00:00" (UTC) but when animating the map based on timestamp the animated time legend displays a range of "03:00 AND 09:00" (which is incorrect).

Is it possible to turn off this auto-correction and display the original time while using Torque?


texts: ArcPy Delete Function not working inside for loop?

texts:
I'm creating a table and then trying to delete that table inside the same for loop as follows:

for fc in fcList:
    arcpy.GenerateNearTable_analysis(inFile, fc, out_table)
    arcpy.JoinField_management(inFile, "Obj_ID", out_table, "In_FID", ["NEAR_DIST"])
    arcpy.Delete_management(out_table)

But for some reason the delete statement isn't working. The table still exists after the script completes 'successfully'. When I just run the delete statement alone afterwards (means outside the for loop), then the table is deleted. But I need to delete it inside the for loop so that I can loop further. I have to put some more statements within the for loop later and the table needs to be deleted first for the upcoming statements to work.

Does the Delete_management statement not work within a for loop for any reason?


texts: Selecting points with MINIMUM distance of 10 km from vector line in country?

texts:
I have a huge amount of points and some lines in my QGIS map.

I want to do a selection just of those points that are minimum 10 km far from the line.

Any idea?

I have no problem doing the opposite (points within 10 km of the line) using buffers - but that's not what i need.


texts: I want to establish connexion QGIS and PostGIS in order to run the pgRouting Layer, how can I do it?


texts: Create random sampling point in each grid's cell

texts:
I would like to know how to create an X number of random points in each grid cell.

I'm using data from a DEM layer, which contains information about elevation. I used this information to obtain slope and aspect. Then I created a grid (each cell with 1000 units size) with the objective of randomly sampling each cell (with one, two points, doesn't matter) so I could obtain an average value of slope and aspect for each location (location = grid cell).


texts: How does OpenLayers determine a layer's projection?

texts:
Could anyone explain the strategy OpenLayers take to determine the layer's projection?

For example, if I didn't explicitly specify any projection for a marker layer, what projection would the layer take?


texts: Converting JSON/GeoJSON files to Shapefiles

texts:
I have a quick question please about converting JSON/GeoJSON files.
I'm trying to convert some JSON files into something I know how to work with (shapefile, R objects or PostGIS). I'm afraid I have no experience with JSON.
I've tried converting them to shapefiles using ogr2ogr, readOGR in R, and FME. I've tried converting them and uploading them to a PostGIS database using FME. And I've tried converting them to R objects using the rjson package. Nothing works.
ogrinfo can't read them and the couple of online JSON/GeoJSON validation tools I've tried say they are not valid JSON files?
Would someone who understands JSON/GeoJSON mind taking a very quick look at one of the files please, and tell me if they're OK? If they're fine, that's great, I can go away and try and work it out. If not, I can save myself a little head scratching. I'm sure there's probably nothing wrong with them.
Here's a link to one of the files in my Google Drive: Opendata Scotland JSON File


texts:
I've recently updated to QGIS 2.14, and have started using the custom variables feature implemented in 2.12. I was wondering whether it was possible to use the custom variable to select a field in the expression builder?

For example, I have a table containing data for three different time periods (AM, IP and PM). If I define a variable "TimePeriod", I'd like to be able to select a field in the expression builder for both labels and style in a way similar to

"@TimePeriod Count"

(where "AM Count", "IP Count" and "PM Count" are three fields in the table)

Is this possible? If so, what would the correct syntax be?


texts:
I’m looking for a set of SVG oil well symbols for QGIS.

Black Bart data: http:// blackbartdata . com / QGIS_Petro_SVG_Symbols

Just directs me to a spam page.

Any other alternatives?

This is a picture of the symbols I ma searching for:


texts: Seeking Free and Open Source GIS Desktop packages?

texts: Can someone give me a list of free and open source GIS desktop packages?


texts: IDE PyCharm settings for developing a QGIS plugin

texts:
How do i setup pycharm to recognize QGIS modules?

IDE: PyCharm
OS: Ubuntu 18.04
QGIS: 3.8

example of imports

this answer is very close to the solution but I cannot add the interpreter location

I've found the location of the plugins folder but I cannot see a python executable here

Here is a screen shot of the "add interpreter" dialogue which I have tried to add the folder location but as you can see it is grayed out. Any suggestions.


texts: Setup PyCharm to develop Python scripts for QGIS?

texts:
I want to set up for the first time an environment with PyCharm to write scripts for QGIS on my Ubuntu desktop and I'm completely lost.

I'm not having problems with the console on QGIS, but I couldn't use PyQGIS on my PyCharm projects the same way I did with pandas package e.g.

Apparently there are some tutorials using Windows but they fail in this case because I really don't know how to adjust for Ubuntu. For example: Set up PyCharm for PyQGIS 3 development. I couldn't find a way to do this with Ubuntu.

Python version: 3.7.4

PyCharm version: 2019.3.3 (Community eddition)

QGIS version: 3.4.10

Ubuntu version: 19.10


texts: split attribute using ArcGIS Field Calculator

texts:
I have a long string in one attribute. I want to be able to extract a small portion of it into a new attribute.
This is an example of the attribute:

"< html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt">   < meta http-equiv="content-type" content="text/html; charset=UTF-8" > < /head > < body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;" > < table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px" > < tr style="text-align:center;font-weight:bold;background:#9CBCE2" > < td >Fly500< /td > < /tr > < tr > < td > < table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px" > < tr > < td >Name< /td > < td >Fly500< /td > < /tr > < tr bgcolor="#D4E4F3" > < td >Notes< /td > < td >< /td > < /tr >  < tr > < td >Source< /td > < td >< /td > < /tr > < tr bgcolor="#D4E4F3" >
< td >Duplicate< /td > < td > < Null >< /td > < /tr > < tr > < td > Type < /td > < td > Wildlife Sensitive Area< /td > < /tr > < tr bgcolor="#D4E4F3" > < td >Start_Date< /td > < td >May 1< /td > < /tr > < tr > < td >End_Date< /td > < td

In bold is what I want to extract. I tried splitting using Right and Left, But it is not always located at the same number of characters from either side. What I would need to do is take everything right of < td > Type < /td > and extract what is between the next  "< td > < /td >" area.


texts: How do I go from a TIFF with symbology to polygons for each label?

texts:
So I am trying to take this TIFF layer and convert it to polygons with their designated label

For example I want a shapefile that has polygons for "-11.374 - 0", "0 - 6", "6 - 13", and "13 to 20". I know that I need to use the tool raster calculator and other Spatial Analyst tools but confused to where to start.


texts: Comparing columns in ArcMap with wildcard statement?

texts:
I need to compare a street centerline name column to a parcel centroid column that holds the street address.

I would like to use a wildcard search so that I can match them by the street name despite prefix/suffix etc.

My goal is to match the corner parcels to the road that it belongs to even though they have multiple frontages.


texts: QGIS street labelling when rotating in composer

texts:
When labeling streets in a QGISmap everything's OK, no label are upside-down.
But I need to rotate at 180° in the composer. Labels of lines don't flip, so they are upside-down!

I tried to "show reversed labels" when a "rotation is done" (sorry my QGIS is in French, words may not be the right ones); and I tried "line orientation dependent position" but nothing happens.

Does someone know if it is possible to flip the lines or the label curved on the line?


texts: Function only mapped over the first image of the collection

texts:
If anyone could help with this please.
I need to map a function over an image collection. The problem is that the function is applied only on the first image of the collection and not on the rest. Here is the code for it.

A = ee.ImageCollection("MODIS/006/MYD09GA"),
    ocean = ee.ImageCollection("MODIS/006/MYDOCGA"),
    geometry1 = /* color: #d63000 */ee.Geometry.Point([10.644168108701706, 64.4188479841206]);

var b1 = A.select('sur_refl_b01')
var oceanb11 = ocean.select('sur_refl_b11')

var filter = ee.Filter.equals({
  leftField: 'system:time_start',
  rightField: 'system:time_start'
});

var innerjoin = ee.Join.inner();

// Inner join applied. The output is an imagecollection.
var innerjoinedMODIS = innerjoin.apply(b1, oceanb11, filter)

var joinedMODIS = innerjoinedMODIS.map(function(feature) {
  return ee.Image.cat(feature.get('primary'), feature.get('secondary'));
})
//print (joinedMODIS)

var dated = joinedMODIS.filterDate('2002-07-05', '2003-07-01')
print(dated,null,'dated')

//PRI for single image
var single = ee.Image (dated.first())
//print (single)

var addPRI = function(image) {
  var pri = single.normalizedDifference(['sur_refl_b11', 'sur_refl_b01']).rename('PRI');
  return single.addBands(pri);
};

//Test the addPRI function on a single image.
var testpri = addPRI(single).select('PRI');
//print (testpri)
//Map.addLayer(testpri,null,'testpri singleimage')

// pri for dated collection
var datedPRI = dated.map(addPRI)
print(datedPRI,null,'datedpri')
print(ui.Chart.image.series(datedPRI, geometry1, ee.Reducer.mean(), 1000))


texts: How to calculate highway lengths using QGIS?

texts:
I need to create a network of highways in 8 states of the USA.

Eventually, I need their length, the speed and if possible average daily traffic volume, I don't know how to extract this information.

I am using QGIS (for the first time today), to create my map.

I already have the map, so my question is how to calculate length of roads?


texts: Selecting points having 1, 3 or more lines intersecting using PostGIS

texts:
Here's the problem : I have two tables, one with roads, the other with points (which are startpoints and endpoints of each road sections, without duplicates).

I want to extract only points that have 1, 3 or more different road sections intersecting.
1 = edge of the network
2 = unecessary nodes between two lines (to erase)
3 and more = network junction

I need those points because of their importance in the network, here's a picture to illustrate the issue :

Basically, I need a PostGIS query that will return point 1 and 3 and leave point 2.

Edit 1 (poor optimization):

DPSSpatial's solution worked fine, even with ST_DWithin. Here's the code :

CREATE TABLE "ReseauRoutier".highways_nodes_clean AS (
SELECT * FROM (SELECT nd.gid, count(
    CASE WHEN ST_DWithin(nd.geom, sg.geom, 0) THEN 1 END) AS count_roads, nd.geom
    FROM "ReseauRoutier".highways_nodes_test nd,
         "ReseauRoutier".highways_segments_test sg
    GROUP BY nd.gid) AS counts
WHERE counts.count_roads > 2);

I think ST_DWithin using a tolerance of 0 may not be appropriate though and might affect the process time on larger datasets...

Edit 2 (also poor optimization):

Just a quick update to show my current code with the ST_Intersects function. Process time is now reduced by half.

CREATE TABLE "ReseauRoutier".highways_nodes_clean AS (
WITH counts AS (
    SELECT nd.gid, count(CASE WHEN ST_Intersects(nd.geom, sg.geom)
        THEN 1 END) AS count_roads, nd.geom
    FROM "ReseauRoutier".highways_nodes nd,
         "ReseauRoutier".highways_segments sg
    GROUP BY nd.gid)
SELECT * FROM counts
WHERE counts.count_roads != 2);

I decided to keep nodes that have only 1 intersecting road as they represent the network's edges.

Edit 3 (great optimization):

I had to edit this post as I think the following is the best solution for this issue. Using ST_Intersects inside the CASE will force PostGIS to test every lines for every nodes, I soon hit the temp file limitation when I used it on large dataset. Thus, my co worker came up with the join idea.

CREATE TABLE nodes_clean AS (
WITH joined AS (SELECT n.id AS node_id, count(w.id), n.geom
    FROM nodes n
    INNER JOIN ways w
    ON ST_Intersects(n.geom, w.geom)
    GROUP BY node_id, n.geom)
SELECT * FROM joined j
WHERE count != 2)

That request takes 1.5 sec instead of 50 sec, which is quite an improvement.


texts: BC30456 'ITableSelection_SelectionChanged' is not a member of 'StandaloneTableClass'

texts:
I am seeing this error BC30456 when building this DLL for VB.Net ArcMap extension. I see the following in the ESRI.ArcGIS.Carto.xml which confuses me with the error.

<member name="M:ESRI.ArcGIS.Carto.StandaloneTableClass.ITableSelection_SelectionChanged">
  <summary>Fires the layer update event. Required when SelectionSet changes.</summary>

Here is the progression

LodeNodes()

cAttachStandAloneTable

MyBase.ITableSelection_SelectionChanged()

Private Sub LoadNodes()
    Try
        lblStatus.Text = "Populating the tree view with the database layers"
        Me.Cursor = Windows.Forms.Cursors.AppStarting

        ' Get relAttachments table and set selection table for it
        Dim relAttach As ESRI.ArcGIS.Carto.StandaloneTable = New cAttachStandAloneTable(mFWS.OpenTable(modDBNaming.GetQualifiedTableName("relAttachments", mFWS)), "FCName")
        mRelAttachSel = relAttach

        ' Load relAttachments into dataset datatable
        PopulateRelAttachmentsTable(relAttach.Table)
        Dim relAttachDV As DataView = New DataView(mDataset.Tables("relAttachments"))

Public Class cAttachStandAloneTable
   Inherits ESRI.ArcGIS.Carto.StandaloneTableClass

   Private _DisplayField As String
   Private _Name As String
   Private _Table As ESRI.ArcGIS.Geodatabase.ITable
   Private _Valid As Boolean
   Private _SelectionSet As ESRI.ArcGIS.Geodatabase.ISelectionSet

Public Overloads Sub ITableSelection_SelectionChanged()
    MyBase.ITableSelection_SelectionChanged()
End Sub


texts: How to remove data points that don't overlap with other data

texts: I am working with a dataset for the whole state, but I want only the data point that overlap with another polygon dataset. How do I get rid of the excess points that don't overlap?


texts: Use GeoServer's REST API on feature names with special characters

texts:
By default, GeoServer comes preconfigured with some example workspaces, datastores and layers. One of these workspaces is called it.geosolutions.

If I try to delete this workspace via cUrl and the REST API:

curl -v -u admin:geoserver -XDELETE  http://localhost:8080/geoserver/rest/workspaces/it.geosolutions?recurse=true

, then cUrl returns

HTTP/1.1 404 Not Found

However, if I run the same command on a workspace that does not contain special characters, such as the pre-configured cite workspace:

curl -v -u admin:geoserver -XDELETE  http://localhost:8080/geoserver/rest/workspaces/cite?recurse=true

, then cUrl returns

HTTP/1.1 200 OK

Hence, I guess that the problem in the former case is due to the . in the workspace's name. How can I properly escape such special characters in the URL?


texts: Are there methods of filtering gps-data with metadata?

texts:
In our project we captured GPS-tracks from a few electric vehicles in the timeframe of a year. In total we captured around 13000000 GPS-points this way. The captured data comes from electric powered cars and we gathered some information from the car as well(power output, average loading time, speed, gear shift, brake state, ...). The Gps-points and the other information is not synchronized. We need a way to filter both the information and the gps so that we can figure out, for example the average output for a part of the road. The data looks something like this:

GPS:

lat1,long1,timepoint1
lat2,long2,timepoint2
lat3,long3,timepoint3

(encoded in NMEA)

Power output:

ev_power1.1,timpoint1.1
ev_power1.2,timpoint1.2
ev_power1.3,timpoint1.3
...
ev_power2.1,timpoint2.1
ev_power2.2,timpoint2.2
ev_power2.3,timpoint2.3
...
ev_power3.1,timpoint3.1
ev_power3.2,timpoint3.2
ev_power3.3,timpoint3.3
...

where timepoints 1.1-1.n are somwehere between timepoint1 and timepoint2.

The final output should look something like this:

lat1',long1',ev_power1',1timepoint1'
lat2',long2',ev_power2',timepoint2'

After the Filtering we want to write it in our Postgis database. We hope to reduce the scattering of our data without losing too much of our extra information. Does anybody have a list of methods I could look into?


texts: Raster DEM interpolation within overlapping vector polygon of water surfaces?

texts:
I am analyzing orthoimage and LiDAR dataset to find out flood areas. I extracted water areas as vector polygon from orthoimage using multiresolution segmentation in eCognition, and generated DTM (1.5 spatial resolution) from LiDAR dataset using moving planes algorithm in OPALS.

Now I am going to extract elevation of water surface areas from DTM and attach to corresponding water features by zonal statistics and also build cross section profile of water surface. One of the drawback of this work is DTM has a lot of big gaps in water areas due to very poor LiDAR returns on water surface.

I am using zonal statistics and cross-section profiles for flood analysis, but getting poor result because of many holes on the water surface in DTM. The profile line was taken from the upper left to the lower right (see the yellow line on the last map).

I tried to interpolate the DTM, but the methods I used interpolate the surface using also the cells of coasts of river resulting false elevation on water surface. I extracted only water areas on DTM and want to interpolate it within overlapping water polygon, so I would get nicer zonal statistics and smoother cross-section profile. Does anybody know the method or tool to do it?


texts: Information about NADADEM files .inc and .inc0 files

texts:
By NASADEM, it is available to download DEM data for worldwide, they sort theirs data into subdirectories hgt_merge, hgt_srtmOnly, inc, inc0 etc.

I am interested in the directory inc and inc0, they contain the tiles as files with .inc extension and .inc0 respectively.

following the documentation:
.inc files are:

  "Subswath radar incidence angle (relative to ellipsoid)"

and .inc files are:

  "Subswath radar incidence angle (local)"

I know already what these angles are, but my question is, what are the shorcuts and the numbers by naming such files 081_080_ss1_d_1_1 etc.
Each tile has 11 own files, for example the tile N52E009 has:

n52e009_081_080_ss1_d_1_1.inc
n52e009_081_080_ss2_d_1_1.inc
n52e009_097_080_ss2_d_1_1.inc
n52e009_097_080_ss3_d_1_1.inc
n52e009_097_080_ss4_d_1_1.inc
n52e009_122_280_ss4_a_1_1.inc
n52e009_138_280_ss1_a_1_1.inc
n52e009_138_280_ss2_a_1_1.inc
n52e009_138_280_ss3_a_1_1.inc
n52e009_138_280_ss4_a_1_1.inc
n52e009_154_280_ss1_a_1_1.inc


texts: How to get a list of features from FME? Python?

texts:
I would like to get a list of geodatabase feature tables from FME.
Is it possible with Python? Or another transformer?


texts:
I have received a request to convert some terrain data to .3arc format. Anybody heard about this format? I did a quick search online but couldn't find anything.
Thanks in advance.


texts: Merging attributes by location in PostGIS

texts:
I am a newbie on PostGIS and I want to transfer an attribute or even a symbology from a layer to another. For example, I have a file with points that each country has a point at the center, I want to transfer the symbology to the country borders' polygon shape file to get a better visualization.
How can I merge attributes or transfer symbology?


texts: Dynamically set default file geodatabase to newly created one in for loop with ArcPy

texts:
I'm trying to write a script for a tool in ArcGIS Pro 2.x for an annual workflow where I create five time-peroid directories (named 1-5) and then I create both a new file geodatabase and an ArcGIS Pro project (using the same name as the folder). I then need to set this geodatabase as the default geodatabase.

My script runs to completion and the ArcPy messages state the default to be the newly created one. However, when I open up each new project in Pro, the newly created geodatabase is not set as the default (nor is it even connected).

# 0. Set-up
import arcpy, os, sys
arcpy.env.overwriteOutput = False

# Read in tool string parameter, year
year = arcpy.GetParameterAsText(0)

# Get the pathname to this script
script_directory = sys.path[0]

# Instantiate ArcGISProject object and assign variable aprx to it
arc_project = os.path.join(script_directory, 'Classification_Project_Template.aprx')
aprx = arcpy.mp.ArcGISProject(arc_project)

# 1. Create new folder for year
new_year_directory = os.path.join(script_directory, year)
if not os.path.exists(new_year_directory):
    os.makedirs(new_year_directory)

# Create 5 time-period directories (named 1 through 5) > within each, make a copy of .aprx, create new geodatabase, and set it as the default
for t in range (1, 6):
    time_directory = os.path.join(new_year_directory, str(t))
    os.makedirs(time_directory)
    aprx_copy_name = str(t) + '.aprx'
    aprx_copy_path = os.path.join(time_directory, aprx_copy_name)
    aprx_copy = aprx.saveACopy(aprx_copy_path)
    # Not sure if aprx_copy will suffice as new project object, but created new one, aprx_new, as ESRI documentation states: "The method [saveACopy] creates a new output project file, but the project variable continues to reference the original ArcGISProject object."
    aprx_new = arcpy.mp.ArcGISProject(aprx_copy_path)
    # Create file geodatabase for each time-period directory
    geodatabase_name = str(t) + '.gdb'
    geodatabase_path = os.path.join(time_directory, geodatabase_name)
    arcpy.AddMessage('Desired default geodatabase_path is: ' + geodatabase_path)
    arcpy.CreateFileGDB_management(out_folder_path = time_directory, out_name = geodatabase_name)
    aprx_new.defaultGeodatabase = geodatabase_path
    arcpy.AddMessage('Default geodatabase is set to: ' + aprx_new.defaultGeodatabase)


texts: 'Select by attribute' using python in ArcMap 10

texts:
I need to select a feature using SQL before zooming to it in python.  What I'm looking for are replacements for the commented lines in this code:

import arcpy
parcel_num = "123456"
mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
# df.selectFeatures("SELECT * FROM LandRecords.GISDBA.Parcels "
#                   "WHERE PARCEL_NUM='%s'" % parcel_num)
df.zoomToSelectedFeatures()
arcpy.RefreshActiveView()


texts:
I'm using R to do some geospatial analysis.

Basically, I am trying to solve for a method of doing poly.counts when my data is too large for my laptop (and R complains about the matrix being too large).

If I have a table I have geocoded, I know I can assign a projection by using

proj4string(my_list_of_things) <- CRS("+init=epsg:26978")

If I wanted to use a for loop in R (gasp!) to pull out a subset (as I have a large amount of data for my puny computer), I am doing something like:

for (i in c(2000) )
{
  print(paste("Postcode: ", i))

  # this subsets things correctly
  a.nam <- paste("subset.",i,sep="")
  assign(a.nam, my_list_of_things[my_list_of_things$postcode == i,])

  #now is the weird bit
  b.nam <- paste("proj4string(subset.", i, sep = "")
  c.nam <- paste(b.nam, ")", sep="")
  # it prints it correctly
  print(c.nam)

  # this bit doesn't work as I expect
  assign(c.nam, CRS("+init=epsg:26978") )

  # I then want to go on and do something like
  nam <- paste("intersections.",i,sep="")
  assign(nam, poly.counts(c.nam, sometable) )
  # to get the number of incidents in the postcode (defined by c.nam) in
  # the overall regions (in sometable)
}

So it prints correctly, but the expectation that within each loop it would run

proj4string(subset.2000) <- CRS("init=epsg:26978")

However it seems to run something equivalent to

'proj4string(subset.2000)' <- CRS("init=epsg:26978")

So instead of 'fixing' the projection for each subset, it creates a variable with a projection defined and the data is "unprojected".

Does that make sense, and if so, how do I solve it?

Based on the example by @SpacedMan, the following example exemplifies my problem:

z = data.frame(w=1:10, x=1:10, y=1:10, z=letters[1:10])
coordinates(z)=~x+y

for (i in c(1) )
{
  message("Postcode: ", i)
  a.nam <- paste("subset.",i,sep="")
  assign(a.nam, z[z$w == i,])

  #now is the weird bit
  b.nam <- paste("proj4string(subset.", i, sep = "")
  c.nam <- paste(b.nam, ")", sep="")
  # it prints it correctly
  print(c.nam)

  # this bit doesn't work as I expect
  assign(c.nam, CRS("+init=epsg:26978") )

}

I'd expect the last "assign" to correct assign the correct projection to, in this instance, subset.1

So even following @SpacedMan's advice and using [i] notation, I have an issue, albeit another one.

for(i in c(1))
{
  test[i] = z[z$w == i]
  proj4string(test[i]) <- CRS("+init=epsg:26978")
}

so in this instance the test object isn't assigned before entering the loop. However I don't know what I should try creating it as, as something like

test <- list()

before the loop throws up an error.

Suggestions?


texts: Generating Contours from csv Concentrations

texts:
I have a shp file showing locations of various wells on a site.  I have linked (joined) a csv file of analytical results to this file.  I am now looking for a method to generate contours based upon the analytical results.  Idea been to identify by contours areas of similar analytical concentrations.
Any ideas?


texts: Correcting ColorInterp on Jpeg2000 NAIP imagery using .aux.xml

texts:
Using GDAL 2.4.0 with OpenJPEG 2000 v2.3.0.

gdalinfo reports all four bands of NAIP imagery as "Unknown" colorinterp. Sample image here: https://prd-tnm.s3.amazonaws.com/StagedProducts/NAIP/id_2015/46115/m_4611501_sw_11_1_20150702_20151207.jp2

This messes up apps like mapserver, which are then unable to serve out color images.

I attempted to add colorinterp information via an .aux.xml sidecar file. Gdalinfo reports that the sidecar file is detected and read, but the colorinterp is still unknown. Sidecar file contents below (what is shown is the entire file).

Googling around, I see reference to the openjpeg2000 driver using the sidecar file for georeferencing, but no mention of using it for colorinterp metadata. Can I fix the colorinterp using a sidecar file?

I know I can fix colorinterp by doing gdal_translate, but I have an awful lot of these files (~4TB), and any gdal_translation produces files 3x - 4x bigger. I want to avoid having the files get that big. I suspect the commercial encoders USDA is using are better than OpenJPEG2000, so I'd love to just leave the files alone and interpret the existing bands correctly. How would I do that?

<PAMDataset>
  <PAMRasterBand band="1">
    <ColorInterp>Red</ColorInterp>
    <Metadata domain="IMAGE_STRUCTURE">
      <MDI key="COMPRESSION">JPEG2000</MDI>
    </Metadata>
  </PAMRasterBand>
  <PAMRasterBand band="2">
    <ColorInterp>Green</ColorInterp>
    <Metadata domain="IMAGE_STRUCTURE">
      <MDI key="COMPRESSION">JPEG2000</MDI>
    </Metadata>
  </PAMRasterBand>
  <PAMRasterBand band="3">
    <ColorInterp>Blue</ColorInterp>
    <Metadata domain="IMAGE_STRUCTURE">
      <MDI key="COMPRESSION">JPEG2000</MDI>
    </Metadata>
  </PAMRasterBand>
  <PAMRasterBand band="4">
    <ColorInterp>Grey</ColorInterp>
    <Metadata domain="IMAGE_STRUCTURE">
      <MDI key="COMPRESSION">JPEG2000</MDI>
    </Metadata>
  </PAMRasterBand>
</PAMDataset>