texts: Why has Web Mercator (auxiliary sphere) become the web map standard?

texts:
I understand what the difference is between the Web Mercator projection and Web Mercator Auxiliary Sphere (WMAS). I also understand that both Google and Esri have adopted this projection as their primary projection for their web-maps which is why we have specialized functions that not only re-project between all projections, but are specific functions for WM, such as webMercatorToGeographic. So ultimately I was wondering why we use the WMAS projection, and the reason it has become a standard in web mapping. Is it a purely a result of two spatial giants moving in that direction or was it just solely because of accuracy reasons?

Additional Links:
Mercator Projection


texts:
I have a big spatial dataset file that I need to import to the QGIS to get the temperature data from each grid. However, the dataset file is in the format of having longitude as first row, latitude as first column. And the temperature data is in the format of a 2D matrix.

Is there a way to read the file to the QGIS? Or is there an easier way to convert the file to QGIS-readable format?


texts: Create fishnet only where cells overlay point feature using ArcGIS for Desktop?

texts: Is there a way to create a fishnet where cells are only generated where they would overlap a point feature in ArcMap 10.3?


texts: Can arcpy.MakeXYEventLayer() cope with *xls files in script?

texts:
I am trying to create a point FeatureClass from an Excel-File using a standalone python-script. I am applying the method "arcpy.MakeXYEventLayer()" which works fine using *.xlsx files. However, I am having problems using *xls files. It also works for *.xls files using the python-window in ArcGIS directly, but not in my standalone script.

For my script I get this error:

  ERROR 000732: XY Table: Dataset C:...\data_for_GIS.xls\sheet1$ does
  not exist or is not supported

Any suggestions?

Here is the code:

import arcpy
import os

arcpy.env.overwriteoutput = True

folder = r"C:/..."
excelfile = "Test.xls"
coordsystem = "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"

arcpy.env.workspace = folder

table = "sheet1"
intable = folder + os.sep + excelfile + os.sep + table + "$"

arcpy.MakeXYEventLayer_management(intable,"GEOLON","GEOLAT","LAYER",coordsystem)


texts: ArcObjects: Alternative workflow to using a SQL query with over 1000 values in IN statement

texts:
I have a ListBox listing the unique values of a column in a feature class; the user can select any number of values and click a button to then zoom to the set of features matching those values. A function builds a WHERE clause that packs the selected values into an IN statement.

This works great when the number of selected values is <= 1000, but if there are over 1000 then no rows are returned. This appears to be an Oracle limitation (the issue does not occur on shapefiles or file geodatabase feature classes), for which some Oracle-specific workarounds are described in this StackOverflow question, but I want to make this function database-agnostic so that it works with shapefiles and file geodatabases as well as SDE geodatabases.

Can anyone suggest any alternative workflows that are 1) implemented in ArcObjects; and 2) database-agnostic; to using an IN statement with 1000+ values? I had thought about creating an in-memory table consisting of the selected values and joining against it but have not attempted this yet, pending any better ideas that may pop up here! If you know of the specific interfaces I should be looking at, that would be very helpful as well.


texts: Aggregate functions in QGIS using spatial intersection

texts:
I am trying to compute a virtual filed using aggregate functions in QGIS (3.6). I want to compute the average on a field in the child layer based on the spatial intersection with the parent layer. Using qgis_sample data for testing and modifying from the example given in the tool I wrote:

aggregate( 'airports', 'mean', 'ELEV', intersects( $geometry, geometry(@parent) ) )

It does not seem to work and only gives NULL back. The goal is to compute the mean on a given field in the child layer (airport) aggregating by spatial intersection with the parent layer (regions), is this allowed/possible with aggregate function?

I've also looked into mean aggregate function by it seems to me that it doesn't accept geometry intersection either for the filter or the group by clause. Or perhaps I am missing the right syntax for it but the documentation is limited.


texts: How to get GetFeatureInfoUrl from multiple layers in OpenLayers 3

texts:
I have multiple layers in a map with layer group now onclick I want to get the GetFeatureInfoUrl for the different layers on which I am clicking but the problem is that I am getting only the Info for the last layer appended not the others.
Here is my code:
 var wms_layer;
    var map = new ol.Map({
        target: 'map',
        layers: [
        new ol.layer.Group({
            'title': 'Base maps',
            layers: [
            new ol.layer.Tile({
                title: 'Water color',
                type: 'base',
                visible: false,
                source: new ol.source.Stamen({
                    layer: 'watercolor'
                })
            }),
            new ol.layer.Tile({
                title: 'OSM',
                type: 'base',
                visible: true,
                source: new ol.source.OSM({
                    attributions: [
                    'All maps © OpenCycleMap',
                    ol.source.OSM.ATTRIBUTION
                    ],
                    url: 'https://{a-c}.tile.thunderforest.com/cycle/{z}/{x}/{y}.png'
                })
            })
            ]
        }),
        new ol.layer.Group({
            title: 'Hisar Data Layers',
            layers: [
            new ol.layer.Tile({
                title: 'Habitations',
                source: new ol.source.TileWMS({
                    url: 'http://192.168.2.2:8080/geoserver/HissarFullData/wms?service',
                    params: {
                        'LAYERS': 'HissarFullData:habitations',
                        'VERSION': '1.1.0',
                        'FORMAT': 'image/png',
                    },
                    serverType: 'geoserver'
                })
            }),
           wms_layer = new ol.layer.Tile({
                title: 'Hisar PWD',
                source: new ol.source.TileWMS({
                    url: 'http://192.168.2.2:8080/geoserver/HissarFullData/wms?service',
                    params: {
                        'LAYERS': 'HissarFullData:pwd',
                        'VERSION': '1.1.0',
                        'FORMAT': 'image/png',
                    },
                    serverType: 'geoserver'
                })
            }),
            new ol.layer.Tile({
                title: 'Hisar Block',
                source: new ol.source.TileWMS({
                    url: 'http://192.168.2.2:8080/geoserver/HissarFullData/wms?service',
                    params: {
                        'LAYERS': 'HissarFullData:block_boundary',
                        'VERSION': '1.1.0',
                        'FORMAT': 'image/png',
                    },
                    serverType: 'geoserver'
                })
            })
            ]
        })
        ],
        view: new ol.View({
            center: ol.proj.transform([75.7217, 29.1492], 'EPSG:4326', 'EPSG:3857'),
            zoom: 10,
            minZoom: 5,
            maxZoom: 17
        })
    });
    var layerSwitcher = new ol.control.LayerSwitcher({
        tipLabel: 'Layers' // Optional label for button
    });
    map.addControl(layerSwitcher);
    var viewProjection = map.getView().getProjection();
    var viewResolution = map.getView().getResolution();
    map.on('click', function (evt) {
        var url = wms_layer.getSource().getGetFeatureInfoUrl(
            evt.coordinate, viewResolution, viewProjection,
            {
                'INFO_FORMAT': 'text/html',
            });
        if (url) {
            document.getElementById('info').innerHTML =
            '';
            $('#myModal').modal('show');
        }
    });

Whenever I click on the map I always getting info of the HissarFullData:block_boundary layer. I am using Layer switcher with this .


texts:
I'm using QGIS 2.6.1 and would like to label US state polygons such that each border is lined with the label of its two states. Something like this:

It appears that the current labeling engine only offers the ability to label above, on, or below the line. This results in something that looks more like this:

What I need is options for interior vs exterior labeling. Is there some way that I can achieve this effect? Any tips?


texts: How to interpolate noise values on a landsat image time series using flag data?

texts:
I need to create a function to interpolate noise values on a landsat image time series by the mean on a time span of 2 (xt+1 and xt-1).

I´m using the fmask product to detect cloud and shadow, then interpolation is applied.

For one time series:

Since c2 is the vector of fmask time series (2 for cloud and 4 for shadow), and t2 the vector of evi time series:

for (i in 2:(length(t2)-1)){ if (c2[i]==2 | c2[i]==4) t2[i]<-mean(c(t2[i-1], t2[i+1]))}

But it is not possible to do this using the calc function of raster package, because it does not works with functions with 2 parameters.

Any suggestion about how deal with this and apply this interpolation for all the pixels of the raster time series?


texts: Retrieving color information about pixels beneath polygon in ArcGIS Desktop?

texts:
I am using ArcGIS for a geology mapping project. I've mapped 500+ features as polygons using a basemap of a planet's surface. I want to get information regarding the "darkness," that is, the greyscale values (0-255) of the pixels contained within the polygons. Ideally I want the descriptive statistics of all the pixels (mean, mode, variance, etc.) contained within the bounds of the polygon.

Does such a tool exist in ArcToolbox, or has anyone ever heard of a plug-in for such a thing?


texts: Removing boundary pixels from numpy array with python

texts:
I need to get rid of the border points in an array which came from a raster. I first converted my GeoTiff to numpy array using:

input = gdal.Open("file.tif")
array = np.array(input.GetRasterBand(1).ReadAsArray())

Then I want to remove the boundary pixels, so that one column and one row (so two of each totaling four lines) are subtracted in the final array. This is because I need to keep applying functions which cannot use these boundary pixels. So I did the following:

window = []
for row in range(1):
    for col in range(1):
        window.append(array[row:(array.shape[0] - 2), col:(array.shape[1] - 2)])

Next I converted this new list into a numpy array via
    newarray = np.array(window)

Originally I had a total of points = 140, with 14 rows and 10 columns. And if I print the newarray to the console my new total is 96 points, which is expected. But if I check the number of rows and columns in this new array via

print newarray.shape[0], newarray.shape[1]

then my resulting structure shows that I only have 1 row and 12 columns. How does this make sense? I should have 12 rows and 8 columns... Is the formula for getting rid of the bounding pixels incorrect, or is there a different and more efficient way to do this and remove the bordering pixels from each side of the array?


texts: How to import and point sample asc files on QGIS

texts:
I am trying to import .asc data into QGIS. These data are produced by the Stage Project in 2000. I want to import it into QGIS to allow me to point sample various climatic variables, especially temperature from it (Link to Stage Three Project data description).

I can’t import the data. The .asc data are a far longer version of what is show below. They are 300 to 600 pages when opened in Word.

sdprecip_cold.ascgcm     v2.0 30k 034
  i  j       JAN       FEB       MAR       APR       MAY       JUN       JUL       AUG       SEP       OCT       NOV       DEC
  1  1     0.032     0.080     0.049     0.078     0.066     0.059
  3  1     0.054     0.080     0.087     0.059     0.043     0.071      0.034     0.048     0.045     0.065     0.063     0.050
  4  1     0.055     0.051     0.045     0.041     0.050     0.068     0.043     0.034     0.083     0.081     0.055     0.058
  5  1     0.049     0.075     0.030     0.062     0.065     0.122
0.180     0.154     0.223     0.209     0.128     0.112     0.091
179 90     0.060     0.104     0.126     0.114     0.116     0.154     0.205     0.247     0.267     0.151     0.123     0.087
180 90     0.074     0.117     0.071     0.133     0.123     0.216     0.165     0.268     0.161     0.144     0.119     0.072

or

grid_warm.asc     v2.0 42k 034
  i  j     grid
  1  1         2
  2  1         2
  3  1         2
  10  1        2
 11  1         2
 12  1         2
178 90         0
179 90         0
180 90         0

I tried importing using Add raster layer and also Translate (Convert format). I then tried removing the many spaces between the cells and also removing the first line of the data but both did not work. I have ensured my driver is set to AAIGrid.


texts: What are last four columns in Royal Mail CSV PAF data used for?

texts:
Royal mail provides detailed address data for the UK in a format called PAF. You can download sample PAF data here.

Most of the field definitions can be found in the official documentation here, but when you download the actual data, there are four extra fields on the end of each line, which are not defined. The screenshot below is from the publically available PAF sample data:

What are these four fields (Thing One, Thing Two, Thing Three, & Thing Four) used for?


texts:
I have downloaded this website template/theme: http://steelthemes.com/demo/html/indusrty-press/

When I look at the code I see this:

They also include gmaps.js

Neither seem to define the colour (or any other colour).  How are they achieving this?

Is there a way to easily download the lat/long for multiple markers from my Google My Maps?


texts: How do I get GeoJSON data from GeoServer into show up on my Leaflet map?

texts:
I'm trying to create a simple "Green" map for my city. The idea being that it is a website with a map and layers you can turn on and off. I am serving the data through GeoServer and am using GeoJOSN to add the data to the map. I checked the URL for the data and GeoServer is serving up the file in GeoJOSN format.

Eventually I would like to be able to let the user to click on the points on each layer to get more information but I need to get the data showing up on the map first. The leaflet code below is pretty simple but I can't seem to get the data to show up on the map. I refereed to the Leaflet documentation but can't figure it out.

I'm guessing I am just leaving something simple out due to my lack of knowledge about Leaflet. Leaflet web mapping ninja's, a little help please.


texts:
I am trying to put together toolbars with model builder tools I created. I followed the instructions on ESRI's page, but it doesn't seem to want to save.

I created the toolbar in ArcCatalog, as well, and it saves fine. ArcMap doesn't recognize that it was created after I close the mxd. I want to have it available for all mxds, and not just the one I save it to.

Edit: The save when I save the MXD, but I want it to save to every MXD I open, not just the one I saved with it loaded.


texts:
I want to use OpenLayers to draw a simple rectangle (based on 4 coordinates) on my map. I would also like to fill the rectangle with a color.

For example (EPSG:4326):

13.40, 52.50
13.50, 52.50
13.50, 52.60
13.40, 52.60

The documentation tells me to use this:

var poly = new OpenLayers.Bounds(0,0,10,10).toGeometry();

How I am supposed to convert this?


texts: How to split linestrings where start/end points don't touch nearby lines

texts:
I am trying to figure out some basic routing with PostgreSQL/pgRouting and I'm struggling a bit.

Our network of linestrings are not topologically correct, tracks that should be connected don't touch, or if they do, it's by chance.

It is my understanding that to make a network topology for routing I will have to split all the tracks at each intersection or crossing. The problem is most the tracks end/start points don't touch, though they are close.

So I'm looking for a way to either split the tracks into segments using a tolerance for points that touch  or to modify the linestrings themselves, moving the end/start point to snap to the nearby line.

But I don't know how to go about that?

I have this query for splitting the tracks, but it relies on points touching. It works for creating segments where the trail crosses another line, perfectly. But the ST_Touches() part is useless when there is no intersection.

CREATE TABLE trail_split_points as
SELECT DISTINCT ST_GeometryN(ST_Intersection(a.track, b.track), 1) as geom
FROM
  trails as a, trails as b
WHERE ST_Touches(a.track, b.track)
OR ST_Crosses(a.track, b.track)
AND a.id != b.id
GROUP BY ST_Intersection(a.track, b.track);

UPDATE

I determined I had no choice but to fix my tracks so they are all topologically correct. It was a bit of a pain to figure out, but I've mostly done that.

The first steep was to get the start and end point of each line and then see if there was a nearby point on another line within <5 meters. If there was then I move that point to that nearby point.

I then use some query to populate a table full of all the points where these lines touch or cross.

Then using that new point table generate new lines for each segment between points. This new table is the one used for routing. The one I created a network map and indexes on for pgRouting.

UPDATE 2

I was snapping the nearby trail to the line, but not creating a new point on that line at the intersection. This is how I managed to do that.

ST_LineMerge(
            ST_Union(
                ST_Line_Substring(track, 0, ST_Line_Locate_Point(track, ST_GeomFromText('".$point['newpoint']."', 4326))),
                ST_Line_Substring(track, ST_Line_Locate_Point(track, ST_GeomFromText('".$point['newpoint']."', 4326)), 1)
        ))


texts: Extracting values of one point shapefile to another point shapefile

texts:
I have a set of points along a river-channel which contain attributes to of the right and left riparian classification respectively. These points almost follow the channel center-line, and have been separated so that I have a point-shapefile with the right's bank classification and another shapefile for the left classification (let's call either shapefile as SHAPEF_CLASS). In the graphic below, the shapefile with the classification (desired attributes) is represented with the bold black dots.

The colored dots represent another shapefile that contains bathymetric data (SHAPEF_BATHY).  I would like to assign the river bank classification attribute from SHAPEF_CLASS to SHAPEF_BATHY. I tried converting the SHAPEF_CLASS to a raster (see below) with the intention to use the Identity tool to extract the classification values. However, the resulting raster does not cover the extent of the river channel so most of the bathymetry point are not covered.

I am using ARCMAP 9.3. I used the default cell size (0.000048), but if I increase the cell size I end up loosing some of the classes from the SHAPEF_CLASS. I thought about converting the SHAPEF_CLASS from point to polyline, and increasing the width of the line so it would cover the entire width of the river channel and then extract the classification. I also have a polygon of the river channel. Not if that would be useful for what I want to do though.

Any suggestions would be welcomed.


texts:
Is it possible to use a Wacom tablet with mapinfo ?

How accurate is it ?

I have a high end machine with dual monitors and I want to digitise field mapping (lots of small polys) quickly.


texts: Raster calculation with multiple NetCDF-files

texts:
I have different climate variables stored in separate NetCDF-files. I like to make calculation e.g. using formulas that consider variable1 and variable2 at each daily timestamp. However, I cannot solve this issue by using the raster calculator although I've aggregated the data to DOY. Is it not possible to use two variables of different datasets at the same time?
I imported the NetCDF as Multidimensional Layers to ArcGIS Pro and aggregated them by recurring days. The new safed format is .crf.
I prefer to solve that issue in ArcGIS Pro.


texts: Projecting points towards selected line layer in PostGIS/PostgreSQL

texts:
I have a point layer (ad.shp) and street layer (st.shp) and I am trying to project the points to the selected street in a distance of 20 meters (snapping all points to the lines within 20 meters) using this approach:

Snapping points in PostGIS

I tried to modify the code given in the mentioned webpage as per points and line layer like this:

-- ad is multipoint layer
-- st is multiline layer
-- Snapping all points within 20 meters
-- Code orginal source: http://blog.cleverelephant.ca/2008/04/snapping-points-in-postgis.html

Select Distinct ON (ad.objectid)
    ad.objectid,
    st.objectid,
        St_AsText(
        St_LineInterpolatePoint(
            st.geom,
            St_LineLocatePoint(st.geom, ad.geom)
        )
    )
From
    (
        Select
            st.geom As Line_geom,
            ad.geom As Address_geom,
            st.gid As Line_id,
            ad.gid As address_id,
            St_Distance(st.geom, ad.geom) As d
        From
            ad,
            st
        Where
            ST_Dwithin(ad.geom, st.geom, 20)
        Order By address_id, d
    ) As subquery;

The error I am getting is "missing From-clause entry for table "ad". Can someone please point to the ambiguity in the above code? I have not tried ST_ClosestPoint() yet.


texts: Problem with double field type in postgis database

texts: I have connection to PostGIS spatial database through ArcSDE, and all double type numeric fields are cutting numbers after comma. Does anyone know how can I fix this problem?


texts: arcpy.MakeFeatureLayer in-memory layer still exists when subsequent step fails during testing

texts:
When I used arcgisscripting to create in-memory output layer with gp.makefeaturelayer, I would put it in a try/except block and delete the gp object in the except piece when the makefeaturelayer would fail.  Now when I am using arcpy, it looks like I need to close the idle window and open again in order to get the in-memory output layer to get deleted.

How could I have the output layer from my makefeaturelayer tool get deleted in a try/except block?  Thanks.


texts: Debugging invalid syntax error in ArcPy script with cursor and where clause?

texts:
I've been trying to write a stand alone script to run any query from python, but I've run into a few problems. At first I tried this:

fc = "U:\RealGDB.gdb\acutalLayerName"
Where = """"TOTAL NUMB FLOORS" = 2""" #TOTAL NUMB FLOORS is the actual field name.
cursor = arcpy.da.SearchCursor(fc, ("'TOTAL NUMB FLOORS'"), Where)

for row in cursor: #python crashes here
    print (row[0])

This crashes python. I've spent about a week trying to fix this and I can safely say that I'm completely lost. Every time I change the sql format to something else - when I add some quotes or brackets - I get an invalid syntax error. What's wrong with my sql; why does it cause python to crash at my for loop?

I imported arcpy and os. The table has about 84.4 records.


texts: Make only particular features on layer selectable

texts:
I'm using a select interaction. I can filter the layers which are selectable using a filter function as the layers option key. I want to also filter which features are selectable.

The documentation says "Which features are selected is determined by the condition option". I am therefore trying to write a function for the condition option which can filter out particular features. But I'm not sure it's possible. I think the quoted documentation may be misleading.

The condition function takes an event argument. I first filter out any events which aren't of the click type, but then I want to filter on features. The event has a target property, but it's not clear to me what this is pointing to. It doesn't appear to be a layer or a feature.

I do have the pixel coordinates and the map coordinates at the clicked point. I could use either of these to loop over features at that point in the layer and, if I deem one of them selectable, return true. But what if there are two overlapping features, one of which I want to be selectable and one not?

I'm guessing this condition option is not what I'm looking for, since it is run once per event, not once per possible feature. Is there a way to achieve what I'm looking for?


texts:
In the description of gdaltransform this example is given:

  gdaltransform -s_srs EPSG:28992 -t_srs EPSG:31370
  177502 311865
  Produces the following output in meters in the "Belge 1972 / Belgian Lambert 72" projection:
  244510.77404604 166154.532871342 -1046.79270555763

If I try this I get

gdaltransform -s_srs EPSG:28992 -t_srs EPSG:31370 177502 311865
244422.238366377 165970.245450285 0

If I write the coordinates into a single point-shape, same result when applying ogr2ogr

ogr2ogr -s_srs EPSG:28992 -t_srs EPSG:31370 out31370.shp in28992.shp

Any ideas what I'm doing wrong?


texts:
I read several threads about it.. How to get $length in meters in QGIS?
Joining attributes from nearest point in QGIS?
QGIS, CRS's and metric measurements of distance matrix and nearest neighbour

and others, here and in other web sites. Nowhere I can find which SRID to use for calculating the distance in meters. My area of interest is the World...anyone please can help me??


texts: Dialogue window for ground to grid tool in arcmap 10

texts:
I am using ground to grid correction tool with some success but in the window I cannot see all the print.

Is there a way to size the window? Does anyone know how to apply ground to grid correction ?
For this is self taught and I see the results but not sure if the workflow I am using is the best approach.

I am using Arcmap 10 on a network.


texts: Options for JS Shapefile to geoJSON conversions that allow multiple input files?

texts:
I am looking for JavaScript options that allow me to take a set of shapefiles (the .shp, .shx, .dbf, and what ever other related files) and convert them into geoJSON.

I am then going to take this geoJSON and save it to the database.

I have found a number of options:

1) https://github.com/wavded/js-shapefile-to-geojson

2) https://www.npmjs.com/package/shapefile

I believe number 2 is my best option but I can't seem to confirm that it will allow me to give it more than one shapefile to convert to geoJSON. Do you know of any libraries that allow for multiple shapefile types to be used as input at one time?


texts: Mosaic multiple multi-band rasters into a single multi-band raster using QGIS

texts:
I have >50 four band rasters that I would like to mosaic into a single four band raster. SAGA mosaic doesn't seem to work since it flattens the rasters out and combines them into a single band mosaic. Building Virtual rasters retains the four bands but does not mosaic so not very helpful across my study area.

Are there any other options for me in QGIS to explore? This seems obvious but I am unable to find something that may work for me.


texts: Toggle ArcGIS Online WebMap Layers using Javascript?

texts:
I am taking an intro course on webmapping and am having some trouble figuring out how to set up a layer list and toggle layers on/off.

I have already preconfigured 2 layers in a WebMap using ArcGIS Online. One is FeatureLayer and the other is a MapImageLayer. My first issue is that my FeatureLayer is not showing up in my legend. The second is that I cannot get the layer list to show up. I would like the ability to turn my layers on/off as needed.

This is what I am trying to accomplish: https://developers.arcgis.com/javascript/latest/sample-code/sandbox/index.html?sample=widgets-layerlist

I've tried copying the list layers function into my code but am struggling to wrap my head around what I need to to get it to work correctly. Also, not sure why my Legend is only showing 1 of my 2 layers.

  <script>
require([
    "esri/views/MapView",
    "esri/widgets/Legend",
    "esri/WebMap",
    "esri/widgets/Locate",
    "esri/widgets/Search",
    "esri/widgets/ScaleBar",
    "esri/widgets/LayerList",
    "esri/layers/FeatureLayer"
  ],
  function(
    MapView, Legend, WebMap, Locate, Search, ScaleBar, FeatureLayer, LayerList
  ) {

    var webmap = new WebMap({
      portalItem: { // autocasts as new PortalItem()
        id: "2ce3aa115bb044a9af0e73929268975e"
      }
    });

    var view = new MapView({
      container: "viewDiv",
      map: webmap
    });

    var locateBtn = new Locate({
    view: view
    });

    // Add the locate widget to the top left corner of the view
    view.ui.add(locateBtn, {
    position: "top-left"
    });

    view.when(function() {
      // get the first layer in the collection of operational layers in the WebMap
      // when the resources in the MapView have loaded.
      var featureLayer = webmap.layers.getItemAt(0);

      var legend = new Legend({
        view: view,
        layerInfos: [{
          layer: featureLayer,
          title: "Legend"
        }]
      });

      // Add widget to the bottom right corner of the view
      view.ui.add(legend, "bottom-left");
    });

    var searchWidget = new Search({
    view: view
    });

      // Add the search widget to the top right corner of the view
    view.ui.add(searchWidget, {
    position: "top-right"
    });

    var scaleBar = new ScaleBar({
  view: view
});
  // Add widget to the bottom left corner of the view
view.ui.add(scaleBar, {
  position: "bottom-right"
  });

view.when(function() {
    var layerList = new LayerList({
      view: view
    });

    // Add widget to the top right corner of the view
    view.ui.add(layerList, "top-right");
  });
 </script>


texts: Relationship between contour interval and spatial resolution of a generated DEM

texts: If I have contour lines with interval say 10, 5, or 2 (etc.), and I generated a DEM from them, what is the relationship between between the contour interval and the resolution of the DEM? Irrespective of the software I am using, whether QGIS, ArcGIS, or any other software.


texts: Trying to build GDAL's ECW plugin in 64 bit Linux

texts:
I'm trying to build the GDAL's ECW plugin in Linux.  I'm using the latest version of the SDK from Intergraph's site (5.0). It fails because it tries to compile against the 32 bit version of the library and not the 64.  Here's the error message where the process dies:

/usr/bin/ld: skipping incompatible
/usr/local/ecwjp2_sdk/lib/x86/release/libNCSEcw.so when searching for -lNCSEcw
/usr/bin/ld: skipping incompatible
/usr/local/ecwjp2_sdk/lib/x86/release/libNCSEcw.a when searching for -lNCSEcw
/usr/bin/ld: cannot find -lNCSEcw

I have pointed ldconfig to look in /usr/local/ecwjp2_sdk/lib/x64/release, but that seems to be ignored.  Is there any way to force it to look in the 64 bit version's directory, or can I get away with changing the symlinks for the 32 bit version and have them point to the 64's instead?


texts: Randomly choosing one from multiple fields in Field Calculator of ArcGIS for Desktop?

texts:
Using the Field Calculator in ArcGIS for Desktop, I want to randomly choose one field from a list of 10 fields (strings).

The Excel equivalent, which is working, is: =INDEX(A2:J2;rnd()*10+1).

How can I do that in VBScript preferably (or in Python)?


texts: OpenLayers select.on event listener fires multiple times (renders multiple)

texts:
I am using hook React. How to render event only once instead multiple times useless?
When I am selecting a feature, it renders multiple times
var select = new ol.interaction.Select();
const translate = new ol.interaction.Translate({
   features: select.getFeatures(),
});

var layer_square = new ol.layer.Vector({
   source: source,
   style: [styles],
});

const layersCollections = [baseSource, layer_square];

const initMap = () => {
   map = new ol.Map({
      interactions: ol.interaction.defaults().extend([select, translate]),
      renderer: 'webgl',
      target: 'map',
      layers: layersCollections,
      view: view,
   });
};

useEffect(() => {
      if (isMounted) {
         map.updateSize()
         select.on('select', (event: any) => {
            var selected = event.selected
            console.log(selected)
            if (selected.length) {
               selected.forEach((feature: ol.Feature) => {
            }
         })
      } else {
         setIsMounted(true);
         initMap()

      }
   },);


texts: ArcMap Dissolve tool only dissolving a few entries

texts:
So I have a shapefile which is divided into 7 habitat categories, with 1845 polygons, and I am trying to dissolve it so I only have 1 line showing each habitat with the total area.

After dissolving, I can see 1804 entries in the attribute table, so it's obviously doing something, but it hasn't dissolved to only one of each habitat.


texts:
This isn't a particularly major issue, but I've had this thing nagging me for a while. When looking at QGIS 3.4's global variables, user_full_name is assigned my given and surname, which includes a diaeresis (¨) somewhere in there. I can't seem to find any way to get QGIS to switch/adapt encoding and display the name correctly. The accented character gets switched to that ubiquitous black rhombus.

I looked for info on this behavior and found moderately similar issues like this and this but it's not exactly the same problem as it's not a geodata encoding issue but rather how QGIS deals with values pulled from the OS.

I have tried unchecking the Override system locale tickbox in Options|General as I prefer to use American English even though my OS is set to Canadian French. No effect to observe.

There's somewhat of a workaround right now which is to fix the name in the Project variable project_author as it is editable (it's filled with the offending global variable by default).

I'm curious to know, though, if someone has encountered this before and if there's a solution in QGIS.


texts: Using Dropbox to Sync File GeoDatabases?

texts:
I've been able to accumulate a large free dropbox account and use it to Sync my GIS databases across several different PC's. I've just started doing this and havent had any problems thus far. But I am getting some red flags of warning that are scaring me.

There seems to be several duplicate files in the GeoDatabases. Usually something like:

a00000017.gdbtable (GIS-HP duplicate Copy)

I can manually delete these files with no problem, but I have no idea why they are popping up.

Also, instead of removing file locks, Dropbox keeps saying it is renaming them. So when I release one lock and gain another, it renames the file instead of deleting it and creating a new one. Although it still does delete the lock file when I have fully disconnected from the GeoDatabase.

Like I said, I haven't had any issues so far, but I am worried about long term data corruption while using Dropbox.

Has anyone used Dropbox and found problems using it to sync GIS data?


texts: I have a great amount of points and I need to show them on map. I am using OpenLayers. Is it possible to simplify point geometry in order to speed up the process of rendering?


texts: QGIS Web Client - postgreSQL Layers wouldn't shown

texts:
since updating my QGIS to version 2.8.3 the postgreSQL Layers wouldn't shown in the QGIS web client.

In the Desktop version the project works very well

I only can see the projecttitle and the groups where the layers in.

Have anybody some idea ?


texts: Avoiding border effect in traffic modeling?

texts:
I'm creating a traffic model (using the Four-Step Travel Demand procedure).

I have a road network (grey lines and grey node) and traffic injector/attractor (orange lines and orange nodes) as it can be seen below:

So the traffic is injected/attracted in my network with those injectors/attractor.

The problem is that I obtain a strong border effect. Indeed I cannot take into account all the trips entering or leaving the network.

The following image show how the traffic flow is concentrated in the middle of my network:

How can I take into account the flow entering and leaving the network ?


texts:
PostGIS is very new to me, and I'm having difficulty performing an operation that probably isn't too difficult. I have vector layers each filled with the outline of a landmass through time. I want to query the distance between these landmasses for every matching timestep. I can do this manually:

select a.gid from_gid, b.gid to_gid, st_shortestline(a.geom,b.geom), st_distance(a.geom,b.geom) as distance
from
(select * from gavdos
where year=3) as a,
(select * from crete
where year=3) as b;

The 'year' value ranges in integers from 0 to 21. It would be convenient to do this query using a FOR loop. I am however highly unsuccessful in doing so, the pgAdmin III SQL thing already presents an error when I say DECLARE i integer;, thus before I am even ready to try and create a loop.

My question is if someone could give an example how I can (1) set-up the loop and (2) how to get the full output (as I suppose I need to update a table at every step in the loop..).

Thanks in advance! :-)


texts:
I've been asked by one of my users if it's possible to find "Rocky Outcrops" in an area. They don't really have a definition for what that means, it's a very open ended question!

My feeling is that this is some sort of "Slope" analysis, combined with a height analysis? We have access to DSM and DTMs from the environment agency as our source data, but this is a new area of GIS for me. Any pointers?

My corporate GIS is Cadcorp, but I can also access QGIS & GRASS if they provide more appropriate functionality.


texts: How to label point features with alphanumeric labels. Only allowing numbers?

texts: New to QGIS. Using point features as markers on a map imported from Google Earth. In my first shapefile layer the alphanumeric labelling worked fine. Added in a second layer and it is only allowing numerical labels? I noticed under layer properties, on the labels tab, the value is shown as 123 Labels (i.e. numeric), whereas for my other layer, the value is shown as abc Labels. Cannnot work out how to correct this?


texts: How do I import/export KML to ArcGIS Desktop?

texts:
I have come across several scripts and tools to import and export KML into ArcGIS Desktop.

Is there a one tool that:

Imports/Exports KML into ArcGIS Desktop at v10
Does not cost money
Ideally maintains attributes between the packages
Ideally does not require additional software to be installed

I would also like to be recommended tools that people have had 1st hand experience with.
ArcScripts, Google, and the blogs I read all talk about a number of methods, but would like to identify a recommended one.


texts: Is it necessary to sanitize GeoJSON/WKT inputs to Leaflet / Terraformer?

texts:
I am currently building a website which includes the following:

a Leaflet map where users can draw polygons using Leaflet-draw (and enter other data using HTML forms)
a paired text field where users can draw polygons on the aforementioned Leaflet map using WKT, which is converted into GeoJSON using terraformer.js
the ability for users to save the drawn polygons/WKT (in GeoJSON) and have them re-drawn later for editing, or...
the ability for users to share a page where anyone can view a map with polygons they have drawn on it (alongside other data).

So essentially: one user can draw a map (using Leaflet-draw or WKT), and many other people can see what that person has drawn, all using Leaflet.
I just added sanitizing of much of the other data using the xss package of node.js, but I am unsure what, if any, sanitization is required for this map. I had thought invalid or unreadable GeoJSON would just produce an empty or default map - but it's possible I am wrong.


texts: STIntersection returning NULL with seemingly valid geography data

texts:
I want to include only the portion of Alaska that's in the western hemisphere. Apparently, some islands are on the other side of the international dateline. I am using STIntersection to intersect the shapefile data I have for Alaska (stored in a geography column) and a text-specified polygon that includes only the western hemisphere. However, STIntersection is returning NULL. The documentation says that it will return NULL if the SRIDs don't match. I checked the SRIDs of both the shape data for Alaska and the text-specified polygon, and they both say 4326, as expected. I'm not sure what to do now. My guess is that I'm not specifying the polygon properly, though it does show up the Spatial Results tab just fine.

This is the SQL to generate the western hemisphere polygon:

SELECT geography::STPolyFromText('POLYGON((-180 0, -1 0, -1 90, -180 0))', 4326) AS s

This appears in the spatial results window as this:

This is the SQL for Alaska. I have a table of shapefile data (from the US census bureau), which I can't easily reproduce here, so you'll have to take my word for it.

SELECT shape
  FROM shapefile_state
 WHERE enum_state = 2

The spatial results window shows this:

If I do a union of those two shapes,

SELECT geography::STPolyFromText('POLYGON((-180 0, -1 0, -1 90, -180 0))', 4326) AS s
UNION ALL
SELECT shape
  FROM shapefile_state
 WHERE enum_state = 2

I get this:

That implies that the polygon should be overlapping Alaska (specifically, the parts in the western hemisphere). And yet the NULL result from the following query:

SELECT geography::STPolyFromText('POLYGON((-180 0, -1 0, -1 90, -180 0))', 4326).STIntersection(shape)
 FROM shapefile_state
 WHERE enum_state = 2

Sanity check:

SELECT shape.STSrid, geography::STPolyFromText('POLYGON((-180 0, -1 0, -1 90, -180 0))', 4326).STSrid
 FROM shapefile_state
 WHERE enum_state = 2

This returns 4326 for both columns.

My knowledge of spatial data is not deep, so I am not sure how to debug this further.

EDIT: here is the shapefile data for Alaska as outputted by STAsText(): http://pastebin.com/AAf5gMSH.


texts: Exporting geodatabase file from ArcMap 10.0 to SQL Server?

texts:
I have a GIS feature class of lines in a file geodatabase that I want to export directly to SQL Server. I only have ArcMap 10.0 to view this file and I can't see how it is possible using that. Usually, I export to shapefile and then upload from there but this is over 2Gb so I want to do it directly.

I already have the database set up in SQL Server, I just want to export directly from ArcMap and create a new table in SQL.

Is this possible?

I have tried using the Feature Class to Feature Class tool, and create a new Database Connection. If I test the connection, it works fine, however when i save the new Database Connection it has a little red cross next to it and still doesn't work. Any ideas?


texts: Suggestions for simple mapping platform for alumni web site

texts:
I'm working on a proof of concept alumni site with a WordPress backend. I want to have a very simple map that shows where alumni are living. I'm thinking something like the Bing API or a simple open source solution - points on a map right now. Not sure yet how I will be getting locations, perhaps a user submission system somehow. Any ideas on good, solid, simple ways to accomplish something like this? I'm on a LAMP stack hosted with WebFaction, and I know just enough PHP to be dangerous.

EDIT: Here's what I ended up with. Jlivni got me looking the Google route and some hacking on this method gave me this:

It's powered on the back end by a publicly viewable Google Spreadsheet. Pretty easy to implement and hack up. And I can pass of the editing of the spreadsheet to just about anyone and they can carry out that task.


texts: Inserting number as well as text in a field in ArcGIS

texts:
I am trying to buffer polygons based on the value of a attribute table, though I need it in centimeter level. There is no option to the buffer tool which uses field to set the centimeter.
Thus, I am trying to to populate the field with value, which is specific buffer distant for each polygon and text "Centimeter" in it, as with this way it seems to work to buffer polygons at centimeter level.
Can anyone share the knowledge how to put number and text in one field?


texts: Counting Files with a Spatial Reference using Python

texts:
I need to display a list of feature classes with an unknown spatial reference, and the number of feature classes with a spatial reference. I think I have the first part right, but having troubles with the second part. I've been playing around with the len() and count() but those always result in an error. The following is part of my code:

for fC in fCList:
     desc = arcpy.Describe(fC)
     spatialRef = desc.spatialReference
     if spatialRef.Name == 'Unknown': # List of feature classes that are unknown
        print "\n" + fC, ("has an 'Unknown' spatial reference\n")
     else: # Number of feature classes that have a spatial reference
        print ("Total number of feature classes with a spatial reference:"), count(fC)


texts: Manually adjusting statistics for raster symbology using ArcGIS Pro

texts:
In ArcMap I can adjust Mean and Std Dev of a raster to get the stretched symbology I want, left side of screenshot.
In ArcGIS Pro Mean and Std Dev are locked.
How can I change them manually?


texts: Setting topology rules to connect short lines aross complex intersections

texts:
I have a polyline SHP with a bunch of very short lines (less than one meter length).

Is there an automated way to force the short lines to connect to a neighboring line at complex intersections?

In the example below, I'd want the highlighted line to connect with the line on the same trend (i.e., northeast to southwest; north is up). Using ArcMap 10.7 Advanced License. I will ultimately need to do this for a few dozen SHP files, and each has thousands of lines in it (so, I'll roll the solution into ArcPy).


texts: OSM: Don't render administrative boundaries

texts:
I have created my own OSM server to pre-generate tiles for later offline use in a desktop app I'm currently coding (I save the tiles in a local SQLite db).
Since I followed blindly a tutorial found online, I'm not sure how to modify what is rendered on those tiles : I would like to have no administrative boundaries shown on the tiles, except country borders.
Right now I have counties, departments, regions etc boundaries and sometimes it's quite hard to differentiate them from paths in the woods.
How can I do that?


texts: Define custom projection in QGIS from some control points


texts: Creating tiled QuickMapServices layer programmatically?

texts:
I’m writing a stand-alone app using Python and QGIS 2.14.9 (upgrading is not an option). Previous to 2.18, tiled services are not supported directly; you need a plugin like QuickMapServices. I have that installed, and it works within the desktop app.

Now I want to use QuickMapServices from my standalone app. I know how to create a QgsRasterLayer, but how can I create a tiled QuickMapServices layer programmatically?

# Create a normal raster layer
layer = QgsRasterLayer(url_with_params, ‘My Layer’, ‘wms’)

# But how to create a tiled QMS layer?
layer = TiledLayer(?)

I can think of two approaches:

Include the plugin path in my pythonpath, and somehow instantiate a TiledLayer directly.
Create a project in the desktop app which includes a QuickMapServices layer, and load that project in my standalone app.

Is there a better way?


texts:
I've been reading through Tyler Mitchel's Web Mapping Illustrated and also downloaded a copy of the itasca sample from maptools.org. Then created some shapefiles in OpenEV and modified the itasca.map file to include those layers and removed all the other layers, but it seems that there's no effect - a blank display
I've also tried to play with the outline color, but nothing. Find the map file below:
find the map file here: http://kompisen.se/xqSjGadQsN

Update:
forgot to specify some details, sorry.

OS: Ubuntu
Mapserver: 5.6.1.1
OpenEV: 1.8 installed via FWTools 2.0.6


texts: How to mosaic multiple modis tiles of several dates using mrt tool?

texts:
How to mosaic multiple tiles of Modis MOD13Q1 in loop. I have about 395 files of each tile from a total of four tiles. I am using linux operating system and therefore use shell script. I need to mosaic the four tiles for each date. I used the mrt mosaic in MRT which automatically took the four tiles and mosaicked it but is there anyway to automate all the (395 X 4) files in a loop?
I carried this out successfully

  mrtmosaic -i mosaicfile.txt -s "1 0 0" -o MOSAIC_temp.hdf

what i need is,

  mrtmosaic -i tempfile.txt -s "1 0 0" -o MOSAIC_$date.hdf ## tempfile.txt temporary text file containing filenames of four files from the same date but from different tiles that are mosaicked and further removed at the end of the loop


texts: Convert coordinates from WGS84 to Gauss-Krüger Zone

texts:
I m currently trying to convert coordinates which come in a csv file in Lon/Lat to Gauss Krüger (Zone 4) using R. My script works perfectly fine, except that the resulting shapefile with the converted points are a little off where they should be (ranging from ~18 to 40m). I guess something is wrong with my spTransform input. As the original CSV didn't come with any metadata, i definded the CRS to WGS84 (the points were recorded with a garmin gps device)

proj4string(results_eDNA)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

where results_eDNA is the dataframe including the coordinates

I then try to transform it with:

coor_eDNA<- spTransform(results_eDNA, CRS("+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000
 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12"))

i think the error comes from the numbers following "+towgs84=", but i have no idea how to figure out which would be the right ones for me. I tried all posibilitys i found on spatialreference.org but with no success.


texts:
Recently I downloaded the new version of eCognition Developer. I remember that before when processing LiDAR data by using algorithms it is pretty easy to generate LiDAR first and last returns.

Before when you chose the algorithm "LiDAR File Converter" there were options under the Converter Parameters such as Convert Parameter and Returns. By setting parameter Returns I could create two new image layers i.e. LiDAR first returns and last returns.

Now with the new version I don't know how to do that because there' no Return parameter listed anymore. I couldn't find any instructions either. How do I generate LiDAR first and last returns in eCognition?


texts: How to install pgrouting 2.x on ubuntu 12.04 LTS

texts:
I try to install pgrouting 2.0 on ubuntu 12.04 LTS with postgresql 9.2

I have ppa in the system by
sudo add-apt-repository ppa:gwibber-daily/ppa

answer:

`You are about to add the following PPA to your system:
 Daily builds of Gwibber trunk
 More info: https://launchpad.net/~gwibber-daily/+archive/ppa
Press [ENTER] to continue or ctrl-c to cancel adding it

gpg: zbiór kluczy ,,/tmp/tmpO_cWbg/secring.gpg'' został utworzony
gpg: zbiór kluczy ,,/tmp/tmpO_cWbg/pubring.gpg'' został utworzony
gpg: zapytanie o klucz 72D340A3 z hkp serwera keyserver.ubuntu.com
gpg: /tmp/tmpO_cWbg/trustdb.gpg: baza zaufania utworzona
gpg: klucz 72D340A3: zaimportowano klucz publiczny ,,Launchpad PPA for gwibber-daily''
gpg: Ogółem przetworzonych kluczy: 1
gpg:               dołączono do zbioru: 1  (RSA: 1)
OK`

next I've typed sudo apt-get update

looks like everything OK, but what next?

on the topic i haven't found much informations.
Thanks for answer.Tom


texts: Enable Add-In Button if One Feature is Selected from Specific Layer

texts:
I am trying to enable an add-in button on my toolbar if only one feature is selected from a specifc Layer. For example, when ArMap loads, the button should be disabled. If the user selects one feature from the AJ_RecommendedImprovementsPoints layer, then the button enables. Here's what I have, but it's not working at all. suggestions??

  Protected Overrides Sub OnUpdate()
    'Get Transportation Improvement Projects
    Dim pFeatLyr As IFeatureLayer
    pFeatLyr = GetLayerByTOC("RV_ImprovementPoints")
    'Enable Button
    If GetFeatSelCount(pFeatLyr) = 1 Then
        Me.Enabled = True
    Else
        Me.Enabled = False
    End If
End Sub

Public Function GetFeatSelCount(ByVal pFeatLyr As IFeatureLayer) As Long
    'Returns # of selected featurues of feature layer as long integer
    'Initialize the required variables
    Dim pMxDoc As IMxDocument
    pMxDoc = My.ArcMap.Document
    Dim pMap As IMap
    pMap = pMxDoc.FocusMap
    Dim pFeatSel As IFeatureSelection
    pFeatSel = pFeatLyr
    Dim pFeatSelectionSet As ISelectionSet
    pFeatSelectionSet = pFeatSel.SelectionSet
    GetFeatSelCount = pFeatSelectionSet.Count
End Function

Public Function GetLayerByTOC(ByVal lyrName As String) As IFeatureLayer
    'This function finds a feature layer based on its TOC Name
    Try
        'Initialize the required variables
        Dim pMxDoc As IMxDocument
        pMxDoc = My.ArcMap.Document
        Dim pMap As IMap
        pMap = pMxDoc.FocusMap
        Dim lyrCntr As Integer
        Dim pFeatLyr As IFeatureLayer
        For lyrCntr = 0 To pMap.LayerCount - 1
            'Ensure that the layer is valid
            If pMap.Layer(lyrCntr).Valid = True Then
                'Ensure that the layer is a feature layer
                If TypeOf pMap.Layer(lyrCntr) Is IFeatureLayer Then
                    pFeatLyr = pMap.Layer(lyrCntr)
                    If UCase(pMap.Layer(lyrCntr).Name) = UCase(lyrName) Then
                        GetLayerByTOC = pFeatLyr
                        Exit Function
                    End If
                End If
            End If
        Next lyrCntr
    Catch ex As Exception
        'Layer not found. Show a message
        MsgBox("Layer " & lyrName & " not found !", vbCritical, "Error")
    End Try
End Function


texts:
I have a geotiff that I've generated the way I'd like it displayed in QGIS. When I load the image as a raster layer, however, it automatically re-scales the pixel brightness values, and the image doesn't display the way I'd like.  When using the GUI, I can click on the layer, then properties. I go to the histogram tab and click calculate histogram, then change the max to 255 (its typically a much smaller value to start, which makes the image look very "washed out"), and the image looks like I want it to..

I'm not sure why it changes it in the first place, but the same thing happens when I add the geotiff as a raster layer via my Python plugin.  However, I can't figure out how to automate this via PyQGIS.. I'm using QGIS version 2.8.6-Wien..

If its important, I make the gtif prior to using QGIS via gdal_translate..

The question is - how can I use PyQGIS to perform this operation so the user doesn't have to do it?  I have PyQGIS automatically adding the raster layer to the project already, its just the forcing it to display the full range of pixel values (within PyQGIS) that I am struggling with right now..

Also - its a 1 band image with a 'singlebandgray' renderer.

Update: the accepted answer pointed to a solution, but for the benefit of others who may be looking for a more direct answer to this particular issue, here is the code as modified to solve my issue:

newRasterLayer = QgsRasterLayer(gtifPathToFname, gtifFnameNoExten)
newRasterLayerProvider = newRasterLayer.dataProvider()

#Force 0 to map to black and 255 to map to white
#Note, this is for a single-band image, so only doing band 0
colorMapList = [QgsColorRampShader.ColorRampItem(0, QColor('#000000')), QgsColorRampShader.ColorRampItem(255, QColor('#FFFFFF'))]

#Jump through some hoops to get the colormap list defined above to
#be used by the new layer's renderer before it is added..
rasterShader = QgsRasterShader()
colorRamp = QgsColorRampShader()
colorRamp.setColorRampItemList(colorMapList)
colorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
rasterShader.setRasterShaderFunction(colorRamp)
pseudoRenderer = QgsSingleBandPseudoColorRenderer(newRasterLayer.dataProvider(), newRasterLayer.type(), rasterShader)
newRasterLayer.setRenderer(pseudoRenderer)


texts: Changing pixel value based on polygon extent in PyQGIS?

texts:
I have a raster layer and a polygon layer. Polygon layer is inside of raster layer and it's extent is smaller than raster layer. Some pixels of raster layer inside the polygon have negative value. I want to change the negative values based on maximum pixel value that is inside the polygon in GDAL or PyQGIS. So In first step I must find the max value based on below code:

r_lyr = QgsRasterLayer('.../aa.tif', 'Raster')
v_lyr= QgsVectorLayer('../bb.shp','bb',"ogr")
v_ext = v_lyr.extent()
#(xmin, xmax, ymin, ymax) = (v_ext.xMinimum(), v_ext.xMaximum(), v_ext.yMinimum(), v_ext.yMaximum())

provider = r_lyr.dataProvider()
r_ext = provider.extent()

rows = r_lyr.height()
cols = r_lyr.width()

block = provider.block(1, v_ext, cols, rows)
lst = []
for i in range(rows):
    for j in range(cols):
        #print (block.value(i,j)),
        lst.append(block.value(i,j))
max_v = max(lst)

In second step it is important to get pixels inside the polygon and then change them to maximum value.

In @xunilk's answer to Changing pixel values of raster data from QGIS python plugin there is a way to change the pixel value but by this answer it is possible to do for full extent of raster layer.

How can I merge my code with @xunilk's code for doing it?


texts: 'No-data' location in predictor variables of Maxent

texts:
In Maxent analysis, is it necessary that all the environmental layers should have 'no-data' at the same places. I have 19 bioclim variables and few others from CGIAR. Both of the sources are not compatible in terms of no-data. The locations of 'no-data' are different in the 2 sources.
Is there any chance of bias, as some of environmental layers have data and others don't at the same location. Does Maxent requires all the layers to have no-data at the same places?


texts: Combining multiple geodatabases into one master geodatabase?

texts:
I currently have geodatabases for each project within a specified area. I would like to combine all of these geodatabases into one master geodatabase. Some of the geodatabases contain some of the same data and if combined into one master geodatabase there would be duplicates. Each database was made with the same template.

How could I go about combining these so the master geodatabase will reflect the updates made to each geodatabase I currently use?


texts: How to update a QGIS plugin for new version (2.2)?

texts:
The updates to QGIS 2.2 have been fantastic and I'm really liking the improvements to the interface, features and performance. However, it seems that some plugins have not survived the upgrade.

cartogram-plugin is one such example, about which there is already a question on this site. My question is not where it has gone, but how I can help get it working again. The latest version of it is, I believe, on GitHub, meaning that a collaborative effort should be able to get it running on QGIS 2.2 is possible.

I am up for giving this a try but would like some advice from the community before I start. Where's the best place to start? There seems to be plenty of advice for people wanting to create new plugins - like the PyPGIS Developer Cookbook and the Plugin Builder plugin. But what about info for people like me who want to upgrade an existing plugin? Sure there is much work to be done in this area, so the answer may be of use beyond cartogram-plugin.

To some extent the question depends on the differences between QGIS 1.x and 2.x. Does that mean this is a good place to start? http://qgis.org/api/2.0/classQgisInterface.html

Many thanks,

Robin


texts: How can I detect overlapping geometries in OpenLayers 6.3.1

texts:
I'm using OL 6.3.1 and have a need to detect when the geometries of polygons overlap or intersect. Specifically I am trying to do something like this when a user draws a polygon:
self.olMap.addInteraction(self.drawInteraction);
self.drawInteraction.on('drawend', function(e){
self.newFeature = e.feature;
var drawnFeatureGeometry = self.newFeature.getGeometry();
var allFeatures = self.displaySource.getFeatures();
for(var i = 0;i < allFeatures.length; i++){
    var thisGeometry = allFeatures[i].getGeometry();
    if(drawnFeatureGeometry intersects with thisGeometry) // pseudoCode here. What can I do?
    {
        self.drawInteraction.abortDrawing();
        self.drawInteraction.source_.clear();
        alert('Features may not overlap. Please try again.');
        break
    }
}

});
I've found examples in previous versions, but those classes seem to have been deprecated.
containsExtent and containsPixel do not meet the use case.
I read the API documentation for Intersects (https://openlayers.org/en/latest/apidoc/module-ol_format_filter_Intersects-Intersects.html) but I couldn't see how it applies.
Any ideas?
Also, abortDrawing() doesn't seem to work as described in the API documentation.


texts:
In an attempt to retrieve all features (point geometry) that fall in a selected polygon, I get the error this.layer.addFeatures is not a function. My code is as shown below:

function createWFSVector(map, WFS_HOST) {
var population = new OpenLayers.Layer.Vector("Population Vector", {
    strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        url: WFS_HOST,
        featureType: "population_layer",
        featureNS: "pe_ns",
        srsName: "EPSG:4326"
    })
});

addMapControls(map, population);
map.addLayers([population]);
}

function addMapControls(map, wfsLayer) {
var select = new OpenLayers.Control.SelectFeature([wfsLayer], {
    box: true,
    multiple: true,
    onSelect: addSelected,
    onUnselect: clearSelected
});
map.addControl(select);
map.addControl(new OpenLayers.Control.EditingToolbar(select));
select.activate();
}

function addSelected(feature) {
selectedFeatures.push(feature);
}

function clearSelected(feature) {
selectedFeatures = [];
}

I can't seem to spot the problem with my code. Any assistance will be appreciated.


texts:
I am current in the process of writing a thesis. Being in the field of GIS it obviously involves some maps ;). While the whole world switched to interactive online maps, Universities still require to hand in a PDF in the end. I am wondering if there are options to integrate a Leaflet map into a PDF. So basically I want to embed a HTML into a PDF.

I also asked a similar question here, but being a more GIS specific on GIS Stack Exchange.


texts: Multipatch multiple slice for creating building floors?

texts:
I'm searching for a method to create multiple horizontal slices at set intervals on a multipatch object.
The slice function in Argis Pro allows only for one horizontal (or vertical) slice to be made, and it is not possible to set a numerical value for the placement of the actual slice.
I'm trying to establish floors in a multipatch object representing a building. I've investigated the "3D Basemap" approach, but this is far too extensive, and not suitable for my project.
There must be a relatively simple way to achieve this either through scripting, or other methods?


texts: QGIS Web Client: Where to upload my .shp files?

texts:
I've installed QGIS Server and QGIS Web Client on an Ubuntu 12.04 Server according to these instructions:
https://github.com/qgis/qgis-web-client

When I load my own project with the Web Client, it shows the correct table of contents but no map, just the box that says "loading map".
A WMS service of my project is running fine and the two sample projects (helloworld, naturalearth) are displayed by the Web Client as well. Looks like the client doesn't find the data for the project.

On the server my folder structure looks like this:

/qgis-web-client/projects/myprojectfile.qgs

/qgis-web-client/projects/shapefiles/myfirstlayer.shp

/qgis-web-client/projects/shapefiles/mysecondlayer.shp

/qgis-web-client/projects/shapefiles/...

The project was created with QGIS 1.8 for Windows using relative paths and this folder structure:

/myproject/myprojectfile.qgs

/myproject/shapefiles/myfirstlayer.shp ...

Where do I have to put the "shapefiles" folder? Do I need to specify that location somewhere?


texts:
I am new to Google Earth Engine and was trying to run a script to get Landsat images of districts of India after having a fusion table of shapes of all districts.
GEE script is as:

var india = ee.FeatureCollection('ft:1UDdgOCf8DoRJ9bVm-UVbR6CqxtkJToLQjTFd0r0Z','geometry')
    .filter(ee.Filter.eq('Name', 'India'))
    .geometry();
var district = ee.FeatureCollection('ft:1_GTWq6xk3hMYl6LelMm8o3VlnoHH0bsFUzLiW3_i','geometry')
    .filter(ee.Filter.eq('district', district_name));

where district_name is a defined variable containing the name of a district

 var india_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
    .filterBounds(india)
    .filterDate('2014-03-01', '2014-09-01')
    .sort('CLOUD_COVER')
    // .limit(500)
    .map(maskClouds)
    .median();
var district_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
    .filterBounds(district)
    .filterDate(year + '-03-01', year + '-05-01')
    .sort('CLOUD_COVER')
    .map(maskClouds)
    .median();

and then after performing few actions having no affect on any of the variables, I am exporting these images to my Google drive.

Export.image.toDrive({
  image: input.clip(district),
  description: district_name + '_summer_' + year,
  folder : 'Project_data',
  maxPixels: 499295920080,
  scale: 30,
  region: district
});

But here after executing and pressing "RUN" button in "Tasks" pane of GEE, I am getting the error "A region must be a GeoJSON Polygon or LinearRing. Got: 'MultiPoint'"

Link to GEE script: GEE script

What could have gone wrong and How do I deal with this?


texts: Adding group layer, renaming it from folder name, then adding layers to it from that folder, using ArcPy/Python?

texts:
I am trying to use os.walk to go through a directory tree, pull out all the .shp files in it, create a grouplayer in the mxd, rename the grouplayer the folder name, then add the layers from that folder to the group. This is because I have a lot of folders with files in them that represent geologic formations and they don't mean anything once seperated from the name of the formation.

From what I read you can't just add a blank group layer in an mxd from python without first creating a blank one and adding that. The problem starts at the for group_name in rename_group_layer loop. It doesn't even make it through to the print group_name test. Basically it just falls out of the loop from what I can tell.

Then I need to figure out how to add the layers from the folder into the renamed group layer. I have two lines commented out that I think are a start, but I haven't gotten it to go through the renaming loop yet. It will rename the first "New Group Layer" if I put that line above the loop to the correct folder name. I am pretty new to python so any help you can give would be greatly appreciated.

import os
import arcpy

# Set workspace environment to a folder
WS = r'c:\test_define\Appalachian Basin'
arcpy.env.workspace = WS

mxd = arcpy.mapping.MapDocument(r"c:\test_define\test_python_add3.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "layers")[0]

print ("#") *80
for (path, dirs, files) in os.walk(WS):
    for name in files:
        if name.endswith(".shp"):
        print path
           print "^^^^^^^^^^^^^^^^^^ path variable"

            try:
                layers = arcpy.mapping.ListLayers(mxd)
                for layer in layers:
                    if layer.isGroupLayer and layer.name == os.path.basename(path):
                        newLayer = arcpy.mapping.Layer(os.path.join(path, name))
                        arcpy.mapping.AddLayerToGroup(df, layer, newLayer, "BOTTOM")
                    else:
                        groupLayer =                       arcpy.mapping.Layer(r"c:\test_define\empty_group_layer.lyr")
                        arcpy.mapping.AddLayer(df, groupLayer, "BOTTOM")

                        layers = arcpy.mapping.ListLayers(mxd)
                        for layer in layers:
                            if layer.isGroupLayer and layer.name == "New Group Layer":
                                layer.name = os.path.basename(path)

                        newLayer = arcpy.mapping.Layer(os.path.join(path, name))
                        arcpy.mapping.AddLayerToGroup(df, layer, newLayer, "BOTTOM")

#mxd.save()

# print messages when the tool runs successfully
            print(arcpy.GetMessages(0))
            #print ("+") *50
        except arcpy.ExecuteError:
            print(arcpy.GetMessages(2))

        except Exception as ex:
            print(ex.args[0])

mxd.saveACopy(r"c:\test_define\test.mxd")


texts: Checking if input raster exist or not in Raster Calculator?

texts:
I am with Raster Calculator in ModelBuilder on ArcGIS Desktop 10.2.

In the Raster Calculator's Map Algebra expression window, I need to check if a Raster dataset has already existed or not.

Can IsNull do this job?


texts: Evaluation and rectification of addesses after geocoding

texts:
I am working on the geocoding of patient addresses for a public health organization.

We use the Normatim (OpenStreetMap) and Google APIs. With the Google service there is no problem so far, but it seems that Normatim, sometimes, has errors in the geolocation of some addresses, particularly with the numbers of the streets.

As you can see in the list, the street numbers are very different (from 2500 to 4300) and should be distributed along the entire street, but Normatim locates them at the same block (image below).

I realize that these are errors in Normatim's data, and I do not intend a solution to this. I am looking for a way to identify these cases (when there is a 'suspicious' grouping of points) to re-process them through Google.

My first approach is to make a heat map, which allows me to visually identify these suspicious clusters and select cases manually.

But  I'm looking for a way to avoid this manual selection, see if there is any way to perform this analysis in a more systematic way. Any ideas or advice?

Note: we must use the Normatim service since it is free and we do not have enough funds (yet) to pay Google (we can only use the 2000 free queries per day that Google offers, but our databases are larger, and growing).


texts: How are Esri's ArcGIS,Google Maps, Mapbox different in their products in providing support maps for Web/App developers?

texts: I am trying to understand their business models. I am not a techie - Angle guy. So kindly explain in layman terms.How is Mapbox innovating?


texts: How to call a ModelBuilder in ArcGIS from a python?

texts:
I want to execute a ModelBuilder from a python script. The ModelBuilder needs two arguments. The model names is m1 and alias mytools.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# man2.py
# Created on: 2017-11-17 09:06:23.00000
#   (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("L:/EPEL/Opérations/MAJ/BgrSigeqRol/4_CompleterTrx/Procedure/Tache9a/AdressesQuébecProductionV3.tbx")

# Local variables:
cf = "C:\\temp\\\database.MDB\\keyfinal"
TEMPO = "C:\\temp\\\database.MDB\\TEMPO"

arcpy.gp.toolbox = "C:\\temp\\\database.MDB\\tool.tbx";

arcpy.gp.m1_mytools(cf, TEMPO)


texts: How to find a feature on map by its field data

texts:
I'm using QGIS 1.8.0 & got much help to find a Featuer using 'Quick Finder' plugin.
Is it possible to locate a Feature (using a similar plugin like'Quick Finder'), on map using its attribute data (Not from Attribute Data Table wizard).

Sandipan Saha.


texts: Retrieving deleted records from feature dataset in enterprise geodatabase?

texts:
I have an enterprise geodatabase with many feature datasets. I was editing one of the feature datasets and accidentally deleted some records from there. I did not realize and stop editing after saving it. When I opened my attribute table, I do not see those selected records.

Is there anyway that I can retrieve those deleted ones in the appropriate data sets?


texts: Get PL/SQL definition of SDE.ST_GEOMETRY functions

texts:
I'd like to learn how to write PL/SQL spatial functions for use in an Oracle ST_GEOMETRY geodatabase. And I figure, the best way to learn something is to start with something that works, thoroughly study it, and then start writing my own functions from there.

So, I'd like to start with an existing SDE.ST_GEOMETRY function. How can I get the PL/SQL definition of SDE.ST_GEOMETRY functions, such as ST_NumPoints (or any other function for that matter)?

I've poked around the web looking for sample SQL for getting the source PL/SQL of functions. There's a lot of talk of packages, etc. and I'm just not really sure where to start.


texts: Stop srcipt working if more than one value are selected

texts:
I have a script tool which has string mulivalue paparameter with a filter list.

For example in this case I want to stop script working after it detect that not only "AAA" parameter is selected but other one or more.
Also it would be useful if script stop after "AAA" is selected in first box and "1" is selected in second.

I could make hundreds of lines like this:

for s in selected:
    if "AAA" in selected and "BBB" in selected:
        arcpy.AddMessage("Don't choose more than one value")
    quit()

But hope that there is an easier way to complete this check.


texts: I want to make it such that my markers at higher zoom levels are replaced by one marker. This is similar to the cluster function, which I have also not been able to make work either. I believe that this would be much easier than the cluster function to make work, and I do not necessarily want to express how many data points there are in any given area. From my understanding the collision function is fairly easy to create but I am unsure how. I'm using custom icons, added from map.addsource(places, {. I've played with the CSS a bit and that has not worked. The collision problems did not seem to exist when I was using style/layer based icons rather than custom ones. My project is adapted from the mapbox tutorial for building stores.


texts: Ohio DNR Oil Gas Lat / Long Conversion in ArcGIS 10.1

texts:
Is anyone familiar with the Ohio Oil Well Search off their website?  For the Lat/Long, they have WH_X and WH_Y.  The numbers in their do not match up with the actual Lats/Longs off the OHio Well Viewer.  This is how it looks.   This is in Belmont County, OHIO (Warren Township)

WH_X    WH_Y
2361000 719070
2359833 714452

Is there a Conversion to the Lats/Longs within ArcGIS 10.1?  Or is there another method within Excel?


texts:
I am trying to do rule-based styling for active drilling rigs in the United States.

Basically, I want to have categorized styling (different colored dots) for all companies (field = "OpCompany") with greater than x active rigs.

For companies with less than x active rigs, I would like them all to be the same color dots.

It is important to note that the data set (picture attached) shows each rig individually, so I am trying to find a way to write an expression/ rule that counts the # of duplicate "Operator Company Name" values in the "Operator Company Name" field.

I have searched extensively online, however cannot find a solution to my problem. Any ideas here?


texts: Update Multiple Columns within One Update Cursor

texts:
Are you able to update several fields within an attribute table at the same item within one Update Cursor command.

Code Example:

with arcpy.da.UpdateCursor (fc, fields) as cursor:
    for row in cursor:
       example2 = (row[0])
    if example == None:
        row[1] = None
        cursor.updateRow(row)
    else:
        variable = dict1[example2]
        row[1] = variable
        cursor.updateRow(row)
for row2 in cursor:
    example2 = (row1[0])
    if example2 == None:
        row1[1] = None
        cursor.updateRow(row)
    else:
        test = somedict[example2]
        row1[1] = test
        cursor.updateRow(row)
 del row

The above code executes but does not update the second column (the first one is updated properly).

The update cursor updates each column based on a dictionary I have created separately. Is this the best way to update a column's value in an attribute table based on a dictionary?


texts: How can I change sort order in layer properties styles tab? The only way I've found is to change their labels and order by it.


texts: Re-licencing for ArcGIS Desktop when no Maintenance purchased?

texts: I have ArcGIS 10.2 in my computer. But recently it suffered OS crash. So obviously I need to re-license the software and I have not subscribed any maintenance package from ESRI. What shall I do?


texts: How to insert data into fusion table using php?

texts:
I am using fusion tables to generate pins over the map. Now I want to insert more coordinates data into the table using PHP but couldn't get any exact idea how to do it?

Can anyone help me for that?


texts: I noticed a change in the Layer Panel of QGIS 3.10/3.12. Vertical grey lines are found at the end of each layer's line. What are they? Do they serve any function? Or are they simply there to help identifying the limits of each layer's line from the right side of the panel?


texts: Caching LiDAR Image Service outside of ArcGIS to send to ArcGIS Server?

texts:
I have very limited programming/coding experience.

We are hosting a number of LiDAR services on ArcGIS Server, and need 24 of these to be cached, from 1:500,000 to 1:4,000.
We initially had published as map services, manually creating cache in ArcGIS after the service was published. However there were blank tiles upon completion (Status read 100%).  We tried publishing the mosaic datasets from the server, from the local desktop, from the NAS... We tried to publish as a cached image service.

We submitted a number of case tickets with Esri and each time they have moved us up to ladder to the next analyst.  The latest response was "No Solution" and suggested to take the open source route to cache our services.  The data itself is fine, and the areas that were covered by the blank tiles cached in ArcGIS are completely visible on the dynamic rendering of the image service.

I need to find a way (outside ArcGIS) to read in a currently published dynamically rendering image service from ArcGIS Server.  I will need to cache this service based on the scale noted above, and save this cache to file location on our NAS to be read in by other applications and users.
I have no idea how arduous a process this will be, nor the requirements to perform such a task.
Apparently it can be done; the image services we host are WMS compatible and support JSON.

I've looked into Tile Mill and MapBox, but honestly wouldn't know where to begin.

If you'd like to look at an example of the image services that need caching, let me know and I can provide a link.

Current version:
ArcGIS Server 10.2.2

We have also tested it on ArcGIS Server 10.3.1 as we are upgrading our machines as soon as we figure this out. But still no solution with the newer version of ArcGIS for Server


texts: Accessing Draw toolbar of ArcGIS JavaScript API without dojo?

texts:
I am building an application using Angular JS, and have been accessing ESRI classes thus far in the following way:

var graphicsLayer = new esri.layers.GraphicsLayer();

I have access to the 'esri' object throughout my application. I am getting this access through my index.html file:

<script src="https://js.arcgis.com/3.12"></script>

Of note, I am not using any dojo code to initialize the map.  Is there any way I can access the Draw toolbar class in the ESRI JS API without having to load it in through dojo?


texts:
I would like to graduate from a column of mine with decimal numbers.
For this I want to apply the Jenks discretisation and exclude zero.
I tried this, but it didn’t work…

Do you have any idea to do that?


texts: I am doing some research and trying to understand how Mapbox compares to GeoServer (or any GIS server software), and whether or not they are used together or if using one precludes using the other. Can you/should you load services stored in GeoServer or similar in a Mapbox map? If you would not use them together, why would you use one vs the other. I'm trying to understand the main differences - it seems to be that Mapbox is more than just an API, perhaps a hosted mapping platform, but how does it interact with GIS servers/databases?


texts:
I'm severely confused by projection systems in GIS softwares and need some help with it.

For my project I need projections which allow me to measure distances and areas in meters or meters squared etc.

From what I know, WGS projections don't allow you to do that, but UTM projections do.

What if I'm using EPSG 32643, WGS 84/ UTM zone 43N? Will this give me measurements in meters?


texts: Adding new field with expression in pyqgis

texts:
I want to use PyQGIS in QGIS 2.18 to add a new field to a vector (line) layer and calculate the bearings for each feature as i have been able to do easily in the field calculator using the expression

concat(floor(degrees(azimuth(start_point($geometry), end_point($geometry)))), '° ',
floor(degrees(azimuth(start_point($geometry), end_point($geometry)))*60 % 60), ''' ',
degrees(azimuth(start_point($geometry), end_point($geometry)))*3600 % 60, '''')

The result of this expression is a string that displays the degrees, minute, seconds value

After a little research, i tried to carry out the task with this script which produced an error.(Actually i was able to use this script to create and update a field that shows the length of lines of the layer, using QVariant.Double and $length for the expression)

from PyQt4.QtCore import QVariant
from qgis.core import QgsField, QgsExpression, QgsFeature

vl = iface.activeLayer()
vl.startEditing()

myField = QgsField( 'Bearing', QVariant.String)
vl.dataProvider().addAttributes([myField])
vl.updateFields()
idx = vl.fieldNameIndex('Bearing')

#next step
exp = QgsExpression('
concat(floor(degrees(azimuth(start_point($geometry), end_point($geometry)))), '° ',
floor(degrees(azimuth(start_point($geometry), end_point($geometry)))*60 % 60), ''' ',
degrees(azimuth(start_point($geometry), end_point($geometry)))*3600 % 60, '''')')
exp.prepare( vl.pendingFields() )

for f in vl.getFeatures():
    f[idx] = exp.evaluate( f )
    vl.updateFeature( f )

vl.commitChanges()

Where am I going wrong?


texts: Summary Statistics tool for Mean Z values

texts:
I'm using ArcGIS 10.2
The summary statistics tool is not functioning in the way that I want. I made an elevation layer that is a DSM, and then in my buildings layer I determined the height for buildings by creating a random points sample. Now that all the random points in each building have an elevation, or z-value, I want to get the mean for each building. I used the summary statistics tool, but it calculates a summarized mean for all 482100 buildings combined. I want the z-mean summarized for each individual building, this is possible because the buildings all have an unique ID number. See Attached Photo,


texts:
I consider myself fairly new to ArcGIS, in that I'm comfortable with all of the basic map-drawing functions. However, there are still things that confuse me.
I am working on a project where I must draw polygons to map landforms over an extensive area. Most features have an ObjectID number. However, there are some that simply have been assigned ObjectID= 0 (zero).
I do not see any immediate issues (ie features are still selectable, editable, etc).
Will there be problems that come up later if I don't address these odd ObjectID's? Why does ArcGIS assign an ObjectID of zero?


texts: Neighbors within distance from every row

texts:
I have a research project of 180,000 entries and I need to calculate the number of neighbors (points from the same table) that exist within a 300 meter radius for each point. My method is very slow, and I'm asking if there is a better way.

My initial thought was to left-join the table to itself, then calculate ST_Distance_Sphere for every row of the result table, then aggregate count the number of entries that are less than 300 meters - but this is turning out to be extremely slow. The initial JOIN in particular is taking more than an hour to complete (it's still running).

SELECT source_id as id, sum(case when dist_meters < 305 THEN 1 ELSE 0 END) as neighbors_within_distance FROM
(SELECT
    a.id as source_id,
    a.biz_lat as source_lat,
    a.biz_lon as source_lon,
    b.id as match_id,
    b.biz_lat,
    b.biz_lon,
    ST_Distance_Sphere(ST_MakePoint(a.biz_lon, a.biz_lat), ST_MakePoint(b.biz_lon,b.biz_lat)) as dist_meters
FROM a LEFT JOIN a as b ON 1=1) as processed group by source_id

I don't currently have any indexes, and the lat/lon columns are stored as floats. I could perform a KNN analysis first to select the nearest 500 points, and then filter the smaller data set down to only neighbors that meet the 300 meter criteria - but I'm not sure how to do that quickly.

Is there a better way to solve this problem?

EDIT : Here are my setup queries and example data rows

DROP TABLE IF EXISTS a;
CREATE TABLE a (
    id              integer,
    name            text,
    street_address  text,
    city            text,
    state           text,
    zip             text,
    county          text,
    telephone       text,
    fax             text,
    email           text,
    web             text,
    type            text,
    num_employees   text,
    sic_1           text,
    sic_2           text,
    headquarters    text,
    revenue         text,
    biz_lat         float,
    biz_lon         float,
    class           text,
    encoding_type   text,
    bbox1           float,
    bbox2           float,
    bbox3           float,
    bbox4           float,
    display         text,
    source          text
);

An example row from the dataset

 ID Company Name    Street Address  City    State/Province  Postal Code County  Telephone Number    Fax Number  Company Email Address   URL/Web Address Company Type    No. of Employees    Primary SIC Code 1  Primary SIC Code 2  Headquarters    Sales/Revenue   lat lon class   type    bbox1   bbox2   bbox3   bbox4   display source
54776   26430. bizname  123 fake    Seguin  Texas   78155   El Paso (830) 555-3910              PRIVATE - PARENT    6   7359    Equipment rental & leasin, nec  Headquarters    1,098,000 (USD) 29.56570181 -97.94911784     place   house  29.56565181 29.56575181 -97.94916784    -97.94906784     123 Fake Seguin Guadalupe County  Texas 78155  United States of America


texts: Toggling layers with arcgis javascript api 3.5

texts:
I am using the two reference samples in the api, toggle and dynamic toggle.  I tried to use the method in the basic toggle but the setvisbiblelayers doesn't seem to like my kml layers and I get a javascript error

Object #<HTMLInputElement> has no method 'setVisibleLayers'

Here is where I add my basemap and kml layers.

            dojo.require("esri.map");
            dojo.require("esri.layers.KMLLayer");

            var layer, map, visible = [];

            function init() {
              var beginExt = new esri.geometry.Extent(-81.26, 40.64, -74.10, 42.29, new esri.SpatialReference({wkid:4326}));
              var revExt = esri.geometry.geographicToWebMercator(beginExt);
              var imageParameters = new esri.layers.ImageParameters();
              imageParameters.layerIds = ['base'];
              imageParameters.layerOption = esri.layers.ImageParameters.LAYER_OPTION_SHOW;

              map = new esri.Map("map",{
                basemap:"topo",
                extent:revExt,
                sliderStyle:"small"
              });

              var kml10url = "https://dl.dropboxusercontent.com/...kmz file";
              var hospital2010kml = new esri.layers.KMLLayer(kml10url, {id:"hospital2010kml"});
              var kml00url = "https://dl.dropboxusercontent.com/...kmz file";
              var hospital2000kml = new esri.layers.KMLLayer(kml00url, {id:"hospital2000kml"});
              map.addLayer(hospital2010kml);
              map.addLayer(hospital2000kml);

This works fine.  I am able to see both of the layers and the basemap.  However when I try to toggle between the layers nothing.  So I tried this, which used to work in 2.6 but no longer works in 3.5.

  function updateLayerVisibility() {
                var inputs = dojo.query(".list_item"), input;

                for (var i=0, il=inputs.length; i<il; i++) {
                    if (inputs[i].checked) {
                        console.log("Input Length " +inputs.length);
                        console.log("Inputs in checked " +inputs[i]);
                        if('hospital2010kml' === inputs[i].id)
                        {
                            console.log('Showing hospital2010kml visibility');
                            hospital2010kml.show();
                        }
                        else if('hospital2000kml' === inputs[i].id)
                        {
                            console.log('Showing hospital2000kml visibility');
                            hospital2000kml.show();
                        }
                    console.log(inputs[i].id);
                    }
                    else{
                        console.log("Input Length Hide " +inputs.length);
                        console.log("Inputs in checked Hide " +inputs[i]);
                        if('hospital2010kml' === inputs[i].id)
                        {
                            console.log('Hiding hospital2010kml visibility');
                            hospital2010kml.hide();
                        }
                        else if('hospital2000kml' === inputs[i].id)
                        {
                            console.log('Hiding hospital2000kml visibility');
                            hospital2000kml.hide();
                        }
                    }
                }

This then gives the error:

  SCRIPT438: Object doesn't support property or method 'hide'

It gives an error for show as well.  Anybody know how to toggle two separate kml layers?  The API says there is a show and hide for kmllayers so I don't know what is up?


texts: How to make a loop for select by attribute?

texts:
I am trying to pull out a sample out of a geodatabase file. I have 100800 sample and I need to select 28 records, basically every 3600, I am going to start with FID=0 and go from there.
Unfortunately I haven't been able to make it work.
this is what I have done so far. anyone can give me a hint, it's gonna be a huge help.

import arcpy, sys, os

Input=sys.argv[1]
field=sys.argv[2]
counts=sys.argv[3]
output=sys.argv[4]

def deleteMe(fc):
    if arcpy.Exists(fc):
        arcpy.Delete_management(fc)
    return fc

cursor =arcpy.da.SearchCursor(Input,field)
for row in cursor:
    a=row[0]
    while row:
        x=a+3600

    arcpy.MakeFeatureLayer_management(Input, deleteMe("lyr"))
    arcpy.SelectLayerByAttribute_management ("lyr", "NEW_SELECTION",a)

#arcpy.SelectLayerByAttribute_management ("lyr", "ADD_TO_SELECTION",x)

arcpy.CopyFeatures_management ("lyr", deleteMe(output))


texts: Displaying bivariate symbology legend in ArcGIS Online as grid

texts:
I'm super excited by the recently added bivariate symbology in ArcGIS Pro. However, when publishing maps and layers, I noticed the legend for my 4x4 grid changes into 16 different lines in the online legend. Not perfect.
Does anyone know of a way to show the legend in ArcGIS Online as it is shown in ArcGIS Pro?


texts:
I have custom infowindows in my map that have embed youtube videos. When I'm working directly in cartodb online, everything loads perfectly.

Within the content section of my infowindow code, I have the following:
< iframe width="640" height="390" src="{{embed_link}}" frameborder="0" allowfullscreen> < /iframe>

I am trying to create a similar map using Cartodb.js (using the api from the Cartodb map I have created with the infowindows that are able to load the video), and the videos were loading fine for a while, but now they have stopped loading. Where the embeded video should be, there's just a blank grey box that says: "The file or directory could not be found."
Is there any reason why the videos might just stop loading when I create a map using Cartodb.js?

Two examples of the embed_link are: "//www.youtube.com/embed/HyCBk9XC2Eo?feature=player_detailpage" and "//www.youtube.com/embed/ZZAFYROf2mY?feature=player_detailpage"


texts: How to pass distances to Multiple Ring Buffer in ModelBuilder?

texts:
I want to use Multiple Ring Buffer in ModelBuilder however I'm not able to pass it 'distances' values. It only accepts one block with a single double value. It's not accepting a list of double values or a string with separated values.

How to pass multiple 'distances' values to Multiple Ring buffer script?

EDIT for PolyGeo:
The model is really trivial:

As long as the distances block is a single double value I'm able to connect it to the script, if it's something else, it doesn't allow me to connect it as Distances. The script itself uses convertValueTableToList function on the entered value. However I can't find an appropriate block with ValueTable type.


texts: Finding percentage of intersection in polygon

texts:
For two given polygons generated in this code, how do I find the percentage of intersection on the rectangle? [The polygon was created using the Alphashape function]

(Linked to Finding points of intersection in polygon? and similar to this question answered: Returning percentage of area of polygon intersecting another polygon using shapely but have not used Shapely to generate polygon for the distribution)

Code implementation:

import sys
from descartes import PolygonPatch
import matplotlib.pyplot as plt
import alphashape
import numpy as np
x = np.random.uniform(41.03641664650488, 59.06406419360736, 1000)
y = np.random.uniform(100.96028140425065, 118.84533708184047, 1000)
test = np.vstack((x,y)).T
alpha_shape = alphashape.alphashape(test, 0.)
fig = plt.figure(1, figsize=(12, 5))
ax = fig.add_subplot(111, aspect='equal')
ax.scatter(*zip(*test))
ax.add_patch(PolygonPatch(alpha_shape, alpha=0.2))
circle = plt.Circle((50,100), 4, color='r', fill=True)
plt.gca().add_patch(circle)
plt.autoscale()
plt.show()


texts:
I found this old post (concerning QGIS 1.x) Is there a Legend Limiter in QGIS?
and was wondering whether this had been solved in QGIS 2.0. Is there a way of only showing the attributes in the legend that are visible on the map?

It seems to me that this is not straightforward in ArcGIS either and I think it would make sense to impement this (in QGIS as I'm trying to use ArcGIS as little as possible). Imagine how useful it would be when you have a geology maps and you don't want to show hundreds of formation but only the few formations you can actually see in the map composer ...


texts: Find the nearest point on a stream network from a given point?

texts:
I have a vector layer for stream network created by r.stream.extract and then vectorized stream file using r.to.vect of grass algorithm in qgis.

I want to find nearest point for a given point on stream network, what  steps should i follow to find nearest point on stream network for this point??

Example below :
Here I need to find the point(X,Y) on stream (line layer) which is nearest to red point.


texts: All 0s result from raster calculator in python

texts:
Raster calculator is not working correctly for a single calculation using a standalone Python script. I am trying to add 1 to a pre-existing raster. It processes, but the resulting tiff file is all 0s. I have used similar code for other rasters and calculations which work, but this one does not.

entries = []
Layer1 = QgsRasterLayer('input.tif', "Lyr1")
print(Layer1.isValid()) # it is valid
ras = QgsRasterCalculatorEntry()
ras.ref = 'lyr1@1'
ras.raster = Layer1
ras.bandNumber = 1
entries.append(ras)

calc = QgsRasterCalculator('"lyr1@1"+1', 'output.tif', 'GTiff', \
Layer1.extent(), Layer1.width(), Layer1.height(), entries)
calc.processCalculation()


texts: Distinct on or unique in QGIS (no postgis)

texts:
I have a shapefile layer which contains mutliple lines (rivers). If i label the layer it will add the river name for each segment of rivers line resulting in mtuliple names along the river. As image

As a quick work round i wanted to use the filter with distinct or unique then use this a label layer. But i run in to a problem when using

select distinct "name"
from layername

The error comes from

The data provider said:
OGR[3] error 1: SELECT from table name failed, no such table/featureclass

How do i find the table name if it's not my layer name?


texts: Modifying existing GRASS GIS tools and running them in QGIS?

texts:
I would like to take an existing tool of the grass tool box and modify the code to the specific needs of my task. As I´m struggling to start some kick-off would be great.

How to access code?

What do I have to consider?

How to access/run the code in QGIS?


texts:
I'm not a GIS specialist at all so this question may seem dumb to you.
I have a GeoTiff base of images (hundreds) I'd like to transform to a WMTS repository to put behind a web server. Web pages are rendered by openlayers.
I'm looking for a way to do GeoTiff to WMTS transformation on linux if possible.


texts: Naming rules for direct connections from ArcGIS to PostgreSQL?

texts:
I don't have problem direct connecting arcgis to postgresql if my database name contains only letters, but when the name contains dashes or underscore e.g:"gisdb-dev", I receive an error when setting up direct connection using ArcCatalog.

  Failed to connect to the specified server. Do you want to continue?
  Invalid SQL Syntax

I'm sure I didn't break any database name restriction on PostgreSQL or ESRI side. I've look around and can't find any solution for this problem.

Had anyone came across this problem and have a solution?

Environment:

Windows 7 Ultimate x64
ArcGIS 10.1 sp1
PostgreSQL 9.1.3 x64
PostGIS 2.0 x64