texts: Using RemapRange for Reclassify in ArcPy gives ERRORS 000622 and 000628?

texts:
I am trying to reclassify a raster, "Conifer_Focl", which lives in inputData.gdb. The max value for "Conifer_Focl" is 33.217064. I'm looking to reclassify as:

I am following the instructions from RemapRange, yet I am receiving the following error:

  File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\geoprocessing_base.py", line
  510, in
      return lambda *args: val(*gp_fixargs(args, True)) arcgisscripting.ExecuteError: ERROR 000622: Failed to execute
  (Reclassify). Parameters are not valid. ERROR 000628: Cannot set input
  into parameter remap.

import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

arcpy.env.workspace = """C:\Users\inputData.gdb"""
inRaster = "Conifer_Focl"
myRemapRange = RemapRange([[0, 1, 100], [1, 2, 85], [2, 3, 75], [3, 4, 65], [4, 7, 40], [7, 10, 20],[10, 33.217064, 0]])

outReclassRR = Reclassify(inRaster, "VALUE", myRemapRange)
outReclassRR.save("reclass.img")

Why does this alternate code work (exported from model)? I don't see an output raster name requested.

arcpy.gp.Reclassify_sa("Conifer_Focl", "VALUE", "0 1 100;1 2 85;2 3 75;3 4 65;4 7 40;7 10 20;10 33.217064 0", "Conifer_Modifier_Model3", "DATA")


texts: Kernel Density Estimation - CrimeStat - lat lon mismatch

texts:
I currently conduct kernel density estimation for a point dataset. For the KDE in addition a reference file is needed, where the coordinates of the lower left and upper right are specified.

The coordinates of my point file are in WGS84, specified with longitude and latitude in the primary file setup. For the reference file, I added the coordinates of the lower left and upper right (longitude for x, latitude for y, as indicated in the documentation).

When I run the routine, I get the following error message:
Irrecoverable error: longtitude/latitude mismatch

The coordinates are specified twice in the same way for the primary and reference file and I already tried different combinations. I also did the method in CrimeStat version 3.3 in order to check bugs of the new version 4.02.

Has anyone a hint?


texts:
I'm using OpenLayers, Geoserver, and PostGIS. For now, everything is fine and happy. I have 2 layers showing on my map through WMS from geoserver. I am taking following steps:

Step 1: inserting data into postgis - done manually.

step 2: if it is a feature in already present layer then reloading the page add that feature on the map, but for adding the table for a layer into geoserver, I am publishing it manually.

step 3: including that layer code into JS manually

I want everything to be done automatically. Since my project is to be handed over to a client and they should be given the possibilitz to add features, layers etc. I think automation is necessary.

So I want:

Step 1: inserting data into postgis - can be achieved by providing a form to user and using a php I can insert them into postgis

Step 2: publishing of layer(selecting a style also) in geoserver - has to be done programmatically

Step 3: inserting this layer code in JS - to be done problematically.

I request your good guidance on this issue. Thank you.


texts: I have a question. I have an AddIn that allows spatial query using map service. I have logic that connects to map service on the client side. I am looking into the possibility of moving this logic to a WCF service so that I could use the same logic for more than one application. I am not sure if I could use rest service from WCF service, based on my search it seems like this is not an option. Is SOAP service is only way to go or if there is a way to use rest service.


texts: Highlighting features that contain data only

texts:
I am trying to highlight individual features in a grid that contain Jobs data from an $.ajax call. I already have a highlight function that highlights the grid cells when they are moused over and I have the cells that do have jobs displaying data for those jobs (in a popup) when moused over.

function highlightFeature(e) {
 var layer = e.target;
 layer.setStyle({
  weight: 0.25,
  fillColor: "#ffff00",
  dashArray: '',
  fillOpacity: 0.5
 });
 info.update(layer.feature.properties);

Now I want to just highlight those cells with jobs data from the get-go so that when the map app is opened, if there are any cells that contain jobs, they will immediately be visible (this will change as jobs are deleted and added by the user). So far, the below code works for adding the popup content for cells with jobs data. How can I use my highlightFeature(e) function to highlight that cell if it's got Jobs data (item.cid === id)? I know there is just something small I am missing in apply the function to that feature (feature.highlightFeature?)

 onEachFeature: function( feature, layer ){
            var id = feature.properties.CID;
            var lat = feature.properties.LAT;
            var lon = feature.properties.LON;
            var popupContent = "";

            JobsArray.forEach(function (item) {
              if(item.cid === id) {
              <?highlightFeature?>
              }


texts: How to display coordinates of a point in a projected coordinate system?

texts:
I want to check the coordinates of a point.

On selecting this point in ArcGIS it displays a latitude and Longitude . In ArcGIS I found an option in 'Identity Tool' that the Latitude and Longitude can be changed into desired units say for eg: Decimal degrees, decimeters etc.

But when the same file is opened in QGIS and the same point is selected it shows me a Latitude and Longitude point but its unit can not be adjusted.

In QGIS, is there a method to change the display units of the point, similar to how ArcGIS does this?


texts: How to refresh the Table of contents in qgis using python

texts:
I have 3 layers in the TOC in QGIS. I want to refresh the layers in the TOC to see what layers are present. Because, when I am executing the following code, I am not getting the layers present in QGIS.

layers = canvas.layers()

So I want to refresh the layers present in the TOC and check. How can I do that in QGIS?


texts:
My mouse no longer allows for panning using a click and hold of the center mouse wheel.  I've tried a new mouse, in case that wheel was broken and no change.  The scroll wheel zooms just fine.

Windows 10. QGIS 2.18.15 or QGIS 3.0.2, neither version works.

Any ideas what settings might have changed?


texts: How do I handle heading in UTM coordinates?

texts:
In a nutshell:
Given a standard (angular) GPS measurement and heading, relative geographic north. If i take my GPS coordinates and turn them into UTM coordinates, do I also have to convert my heading?
I can't find this topic brought up anywhere, so I'm guessing that Im missing something important here.

My reasoning:
As the UTM coordinates deviates the further away you get from the center planar projection point (centered around the prime meridian), the heading must also deviate. Since the prime meridian is the only line that is parallell to the geographic north (except on the equator), the error in heading would increase the further i move away from the prime meridian.

Look at the above image, showing a very large plane projected on a sphere. The gray and red latitudinal lines represents geographic north in point p1 and p2 respectively. On the prime meridian, in p1, the geographic north reading is the same in angular and UTM coorinates. When moving away from the prime meridian, as in p2, the UTM y-axis does not line up with my GPS reading of geographical north. The error is thus the angle between my UTM y-axis and my geographic north.

My solution:
Look at the second image. Subtending angle c would be 90-latitude, subtending angle a would be longitude - UTM prime meridian. Spherical angle B is 90 degrees (orthogonal spherical coordinate axis)

From the cosine rule it can be derived that

cos(a)*cos(B) = cot(c)*sin(a) - cot(C)*sin(B)

Thus, spherical angle C becomes

C = atan(sin(B) / (cot(c)*sin(a) - cos(a)*cos(B)))

And my final error in heading becomes 90-C degrees.

Am I correct in that there is an error present when looking at geographical north from a planar projection, and would this be it?


texts:
I'm looking for inspiration and possible tools on how to display cross-sections on a (dynamic) webpage. More specifically I have to show on these cross-sections how underground pipelines are located  compared to the surface. Apart from the 2D information on the surface (high resolution) we will also have access to lidar data of the surface around these pipelines.

I'm basically considering two paths. One is to make a true cross-section, like eg this example:

Another option would be to create a small strip of land, a bit like this image (but with a much more boring look, but and showing the depth of the pipelines also. . This has the advantage that the surrounding terrain is also shown, which is useful as the position of the pipelines is not always very precise.

Although I have some ideas on how to create the first one (basically making a projection along the transect), I wonder if any tools exist that could create images like the second one. And in general I'd be interested to know if anyone has already knows of (or made) similar applications.


texts:
I am trying to use a ESRI basemap satellite image in order to obtain forest acreage values from a tax map parcel (image included). I have a parcel that contains farmland and forested areas, and I want to find out how many acres of forested land that specific parcel has. I tried the measurement tool in ArcMap but I am unsure of how accurate that is.

Is there any other way to do this?

[


texts: Is there a piecewise-UTM coordinate system?

texts:
Is there an EPSG code, WKT or a proj4 string, for a single global piecewise-UTM coordinate system such as MGRS (the NATO military grid reference system)?

I suppose this would work similar to the HEALpix projection?


texts: PostgreSQL/PostGIS 12. shp2pgsql.exe only blinks

texts:
I am running PostgreSQL/PostGIS version 12 (win10/64 bit) attempting to import shapefiles by means of shp2pgsql.exe.
Really nothing happens. The GUI only blinks a fraction of  a second. Then nothing.
What am I doing wrong?


texts: ArcGIS splitting tool does not update lines properties

texts:
I wanted to split some lines I drew in ArcMap using the "Split" tool during editing.

What I was expecting is that once I split the line in, e.g. 10 segments, for each one of them ArcGIS updates Feature Attributes such as the objectID and the segment length, but I noticed these parameters were not modified.

Why does this not occur? I though the object ID was supposed to be unique for each feature in a layer or shapefile, and that if I splitted a line, the attribute "Shape_leng" (whose name was given by ArcMap itself, I did not created this feature) of each of them would have been updated (continuing with the previous example, I was expecting segments whose length is one tenth of the original one)


texts: Adding ETag to response headers of Geowebcache

texts:
Serving tiles from Geowebcache, I need to add ETag parameter to response headers. How can it would be happen?
Cache-Control: max-age=2629743, must-revalidate
Content-Length: 541947
Content-Type: application/vnd.mapbox-vector-tile
Date: Sun, 16 Aug 2020 07:16:14 GMT
Expires: Tue, 15 Sep 2020 17:45:17 GMT
geowebcache-cache-result: MISS
geowebcache-crs: EPSG:900913
geowebcache-gridset: EPSG:900913
geowebcache-tile-bounds: -20037508.34,-20037508.34,20037508.34,20037508.34
geowebcache-tile-index: [0, 0, 0]
Last-Modified: Wed, 15 Jul 2020 09:31:38 GMT
X-Content-Type-Options: nosniff
X-Frame-Options: DENY
X-XSS-Protection: 1; mode=block


texts:
I'm using FME Desktop and Server 2014.
At the moment my workspace looks like the following. It reads a .csv file and transforms this to a ESRi Shapefile. This is working. I can use the shapefile with QGIS or ArcGIS Desktop.

The next step, which I want to perform:
Create a WMS Service Link, which uses my FME or ArcGIS Server.
Because FME provides the functionality to publish the workspace as an OGC Web map Service, I would assume that this works.

I've already tried to publish the workspace as shown in the first pic. The WMS Service is up and running and I can load the WMS Service to ArcMap, I can query the WMS link and the layer is shown. Unfortunately when I want to add the layer to the map, I get the error "A selected item could not be added to the map. Could not connect to server".

So my first question: Is there any error with my workspace?
Second "question": If searched for some solutions and I found that I need to have a result as PNG, JPEG, TIF.
So my question would be: What do I need to add/edit in my workspace to get my idea to work?


texts: ArcGIS 10.3.1 - Cutting 2D sections from 3D Polygons

texts: I have a bunch of 3D Multipatch features that represent the internal rooms in a large building.  I need to be able to cut 2D sections out for any given elevation, ie make a plan view for any given elevation.  I made a Multipatch cutting plane at one of the desired elevation but looking through the 3D Analyst tools I'm at a loss as what I could do next.  It seems it would be simple, any ideas?


texts: Reading and writing existing KML files using simplekml python module?

texts:
I am trying to add an image icon overlay into an existing KML file (that contains polygon layer from shapefile) without having to overwrite the file. I am currently using the "simplekml" python module which is giving me what I want.

The following code:

import simplekml

kml = simplekml.Kml()
screen = kml.newscreenoverlay(name='ScreenOverlay')
screen.icon.href='http://simplekml.googlecode.com/hg/samples/resources/simplekml-logo.png'
screen.overlayxy = simplekml.OverlayXY(x=0,y=1,xunits=simplekml.Units.fraction,
                                       yunits=simplekml.Units.fraction)
screen.screenxy = simplekml.ScreenXY(x=15,y=15,xunits=simplekml.Units.pixel,
                                     yunits=simplekml.Units.insetpixels)
screen.size.x = -1
screen.size.y = -1
screen.size.xunits = simplekml.Units.fraction
screen.size.yunits = simplekml.Units.fraction
kml.save("ScreenOverlay.kml")

Is there a way I can read and write an existing KML file without having to recreate a new one. Is there a function that allows me to modify and save the KML file?


texts:
I am trying to port some code from an older version of QGIS. I have managed to get the map setup fairly to my liking, and can export some content via a command like

iface.mapCanvas().saveAsImage("example.png")

However, I want to export to PDF, and also add a legend for my graduated color scale.  I've been rummaging around the documentation online for some time but I can't figure out how these new classes work together.

@JoshC helped me to understand that I need to use the QgsLayout and QgsLayoutExporter classes, not mapCanvas, however I'm not clear on how to use those classes to get out a PDF with a legend.


texts: I have two shapefiles, Country and Blocks.  I need to run a very basic selection to find out what blocks intersect the country shapes.  But, the country shapefile has 250+ countries and I need to do this for each row.  The output is a list (ideally a .txt file) of one of the attributes of the block shapefile.  I'm trying to automate this so I don't have to query each country then run the intersection manually.  I'd like to run a Python script or make a model builder, but I don't know how to automate going from one row to the next.


texts: Get 'layers' name form GeoServer link and assign it as a variable

texts:
Using Leaflet and GeoServer.

Basically, I want to be able to add a new layer with a different bounding box that is user defined by selecting an area on the map. I want it to work with all layers available... so if I could just get the 'layers' name and update it when a new overlay is turned on/off

User turns on a tileLayer with many points on the map. User selects an area on the map with points. The selection area is converted to new bounding box dimensions. A new layer of the same map is added but with with the new bounding box.

I have everything except how to get the name of the current tile layer so that I can make the GeoServer call. It would look like this:

var currentLayerName;

  function addGeoServerLayer(){
       geoServerLayer = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: **currentLayerName**, **<<<<<<<<<<<<<<<  I NEED THIS**
            format: 'image/png',
            transparent: false,
           CQL_FILTER : 'BBOX(geom, ' + boundingBox + ')'
        }).addTo(layersGroup);

    }

I want to be able to pull the 'layers' field from this list of tileLayers in my project (leaflet project). this way, after the user selects an area, I can make a geoserver and return the same map with new bounding box values and I can do it no matter what layer is selected.

        var layer1 = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: 'test:layer1',
            format: 'image/png',
            transparent: true
        });

        var layer2 = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: 'test:layer2',
            format: 'image/png',
            transparent: true,
        });

        var layer3 = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: 'test:layer3',
            format: 'image/png',
            transparent: true
        });

        var layer4 = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: 'test:layer4',
            format: 'image/png',
            transparent: true
        });

        var layer5 = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: 'test:layer5',
            format: 'image/png',
            transparent: true
        });

        var layer6 = L.tileLayer.wms("http://12.235.22.99:8080/geoserver/ows?", {
            layers: 'test:layer6',
            format: 'image/png',
            transparent: true
        });


texts: Can SRTM preserve topological feature well?

texts:
There are studies on SRTM data vertical elevation accuracy, which seem to indicate that there are around an average of 5m vertical error when compared with a DEM derived from cartographic maps.

Obviously, with that kind of error, one can't rely on SRTM to do accurate cut fill calculation. But SRTM can still be useful for preliminary site study, as long as the relative contour lines and essential topological features still remain intact.

My site is  an area covered with forest, nearby to a small town ( population around 2000), somewhere in the middle Borneo. I will need to do preliminary site study, like prelim. proposal of roads, drains and important water detention facilities like (man-made) ponds.

Is SRTM good enough in preserving the relative contour lines and essential topological features?


texts: Finding source code of road graph plugin in QGIS?

texts:
I'm working on a project that I have to model a highway map into a network and apply some network analysis algorithms in that network to get insight into.

I found that QGIS Road Graph Plugin routing pretty well only with raw map data. So I searched for source code of this plugin but it seem that the SVN repository hasn't available any more.
I searched in QGIS Github repository, but still cannot find.

Do you have any idea about where they put that code in?


texts: Error in running WhiteboxTools in command line

texts:
I am trying to perform IDW interpolation in WhiteboxTools in my LAS files to create raster image through command line:
C:\WhiteboxTools_win_amd64\WBT>whitebox_tools -r=LidarIdwInterpolation -v ^
More? -wd = "E:/Processed_lidar/norm_thinned_01m/norm_thinned_01m/" -i=GO72A-3.las -o=idw.tif ^
More? --resolution=0.1 --weight= 2.0 --radius=0.3"

And getting following error:
thread 'main' panicked at 'Error parsing -weight: ParseFloatError { kind: Empty }', src\tools\lidar_analysis\lidar_idw_interpolation.rs:326:26
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

I am using windows 64 bit OS.


texts: Changing marker shape to circles in Leaflet

texts:
I was wondering how to change the default marker shape in a Leaflet map to circles

$.ajax("data/StatesPopulation.geojson", {
        dataType: "json",
        success: function(response){
            L.geoJson(response, {
                onEachFeature: onEachFeature
            }).addTo(mymap);
        }
 });


texts:
I have Postgre database with network table that contains country, street name, numbers from left – to left, numbers from right-to right.

I finished geocoding function to use like this: user fills 3  text boxes for country, street and number  and then I have query like this:
Select * from network where country like (text box country) and street name like (text box street) …etc …
For house number I have interpolation and the geocode function returns lon/lat.

Everything is fine with this function.

But I want to do geocoding like Google, user has only one text box and user can insert whatever he wants, and Google knows what to search.

For instance if user insert “Via Romma Palermo” or “Via Roma Palermo” or “Via Roma 16 palermo” Google always has the right answer and is very fast.

What is the trick? How to do fast search, do I need some algorithm for this or what? Do I need separate table for all addresses? How to handle errors in typing Romma, Roma …etc?
How to know what is first country, street, what field in the database to search for?


texts:
Using the following code, how do I get the distance of the route directions?  I want to run some math against the length.

      container: 'map', // container id
      style: 'mapbox://styles/mapbox/streets-v9',
      center: [-74.50, 40], // starting position
      zoom: 9 // starting zoom
  });
  map.addControl(new mapboxgl.NavigationControl());

  var directions = new MapboxDirections({
    accessToken: mapboxgl.accessToken,
  })

  map.on('load', () => {
    directions.setOrigin([locations[0].longitude, locations[0].latitude]);
    directions.setDestination(userLocation);
  });

  map.addControl(directions, 'bottom-left');```


texts: GEE: Getting empty columns when mapping over imageCollection with daily images to reduceRegions

texts:
I'm trying to obtain daily NO2 pollution values for global district areas. I fail to get any of the bands of my images returned as columns. When I reduceRegions over only one image from the collection, I manage to obtain the bands of the image as columns.

Here's the code:

var districts = ee.FeatureCollection(gadm);

var no2Days = ee.ImageCollection("COPERNICUS/S5P/NRTI/L3_NO2")
  .select('NO2_column_number_density')
  .filterDate('2020-04-26', '2020-04-28');

var scaleNo2 = no2Days.first().projection().nominalScale();
print(no2Days.first());

var no2FirstDay = no2Days.first();
print(no2FirstDay);

var no2FirstDayRed = no2FirstDay.reduceRegions({
  reducer: ee.Reducer.mean(),
  collection: districts,
  scale: scaleNo2,
 });
print(no2FirstDayRed.limit(5));

var dailyAverage = no2Days.map(function(image){
  var time = ee.Date(image.get('system:time_start')).format('YYYY-MM-dd');
  var districtMeans = image.reduceRegions({
    'reducer': ee.Reducer.mean(),
    'scale': scaleNo2,
    'collection': districts,
  }).map(function(feat){
    return ee.Feature(feat).set('Date', time);
  });
  return districtMeans;
}).flatten();
print(dailyAverage.limit(5));

Export.table.toDrive({
  collection: no2FirstDayRed.limit(5),
  description:'no2Test',
  fileFormat: 'csv'
});

The script can be found here: https://code.earthengine.google.com/5fd4bda5eb3db4ebd552fc2c66439edc


texts: How to refresh the fields in a popup for Operations Dashboard?

texts:
I have a hosted feature layer that I had added to an Operations Dashboard.

Then I removed 3 fields.

Within the map viewer, when I click on a point, the popup displays correctly with the 3 fields removed.

Within the operations dashboard map, when I click on a point, the popup still displays the old fields.

How do I refresh the Operation Dashboard with the changes?


texts: How can i change vertex option programmatically using python in QGis?

texts:
I want to activate the option "Show markers only for selected features" in Setting-> options->vertex markers ,using python.
does anyone know if i can do it programmatically ?
i attached a photo of the ui and i market the option,


texts: Upload raster data using ArcGIS API for Python?

texts:
I have been using the code below to import shapefiles to ArcGIS Online. It was taken from an ArcGIS tutorial linked here: https://developers.arcgis.com/labs/python/import-data/
from arcgis.gis import GIS

# Log into your ArcGIS Online account
gis = GIS("https://www.arcgis.com", username="your_username", password="your_password")

# Make a dictionary to store metadata for the shapefile
shp_properties = {
    'title': 'Los Angeles Hub Datasets',
    'tags': 'testing, this, with, python, code',
    'type': 'Shapefile'
}

# Load the data as a feature layer with a ContentManager by calling gis.content.add which returns an Item object.
data_file_location = '/path/to/Parks_and_Open_Space.zip'
shp_file = gis.content.add(shp_properties, data=data_file_location)

# Call the publish method to publihs the Shapefile, which returns another arcgis.gis.Item instance for the feature layer.
shp_feature_layer_item = shp_file.publish()

# Identify the URL
print(shp_feature_layer_item.url)

But the tutorial states it can only import CSVs, XLS,  GPX, GeoJSON, and shapefiles. Is it not possible to use Python to upload raster data to ArcGIS Online?


texts:
I have point layer and a network (segments) stored in a geopackage and using QGIS. The number of points and their ids are always different.

What I need is to get a route from first to last points but through several stops in a specific order.
I am thinking in a script or algorithm where to select an order column and then do the shortest path between first and last points setting the rest of it as stops ordered and get a single layer with the entire path.
I've tried the analisys tools of QGIS and plugins but I can't find a tool with this functionality.
I have also tried the v.net.salesman that computes a route through all the points but it isn't what I need, It can't assign an order to the stops and  the route should not be cycle.
Is there any way, using QGIS or PyQgis, to do that?


texts: Assign different colours to intersection between two polygons and difference area

texts:
I am quite new to GIS using R.

I have two sets of data, upon which I could create two polygons like this:

What I want to do is colour the common area (inner 'circle') and the space between borders. However, if I fill the two polygons directly, the colour for the outer polygon simply overlayed the inner polygon. Is there any way to do this nicely?
My codes are:

ggplot()+ geom_polygon(data = inter_Coords, aes(x = inter_Coords[,1], y
= inter_Coords[,2],alpha = 0.25)) +
geom_polygon(data = concave_pts, aes(x = concave_pts[,1], y =
concave_pts[,2], alpha = 0.25)) +
geom_polygon(data = diff_Coords, aes(x = diff_Coords[,1], y =
diff_Coords[,2], fill = 'yellow', alpha = 0.25))

Results are:


texts: cluster features styling (by value return of the function)

texts:
Whether it's possible to insert the return value of the function to the field of property?
I would like cluster features styling depends of average temperatures. My code:

    var highRule = new OpenLayers.Rule({
filter: new OpenLayers.Filter.Comparison({
    type: OpenLayers.Filter.Comparison.GREATER_THAN,
    property: "${val}", <---------THIS------------------------------------------>
    value: 15
}),
symbolizer: {
    fillColor: colors.high,
    fillOpacity: 0.9,
    strokeColor: colors.s_high,
    strokeOpacity: 0.2,
    strokeWidth: 15,
    pointRadius: "${point_radius}",
    label: "${val}", //°C",
    labelAlign: "cm",
    labelOutlineWidth: 0.5,
    fontColor: "#00376d",
    fontOpacity: 1,
    fontSize: "10px",
    fontWeight: "bold"
}

    });

    var style_temp = new OpenLayers.Style(null, {
rules: [singleRule, lowRule, middleRule, highRule],
context: {
            val: function(feature) {
                if(feature.cluster.length > 1) {
                    var sumValue = 0;
                    for (var i=0; i<feature.cluster.length;i++)
                    {
                        sumValue += parseFloat(feature.cluster[i].attributes.tzewn);
                    }
                    var averageValue = (sumValue/feature.cluster.length).toFixed(1);
                    return averageValue+"\n°C  ";
                                                }
                else
                    var Value = feature.cluster[0].attributes.tzewn;
                    return Value+"\n°C  ";
                                    },
            point_radius: function(feature) {
                if(feature.cluster.length > 1)
                    return (feature.cluster.length/9) + 10;
                else  return  8;
            }
    }}  );

This code doesn't work, so syntax ---> property: "${val}" <--- is not correct.


texts: Conditional combining of multiple rasters to produce one output using Python?

texts:
I have two rasters for same spatial extent and need to compare both rasters cell by cell and produce one output raster. Each has 10 land cover classes. Based on those 10 classes the comparison should be  100 conditional statements.

Here I wrote two statements and did a test run. It gives me only a output based on the last conditional statement. Which is output raster has cell values only 21.

outCon = Con (((inRaster1 == 11) & (inRaster2 == 11)), inRaster1 + inRaster2, 11)
outCon = Con (((inRaster1 == 21) & (inRaster2 == 11)), inRaster1 + inRaster2, 21)
outCon.save ("G:\\LC_Transitions\\con_test")

Any suggestions on what I did wrong in here?

Also, if I want to run a statement for multiple conditions together (for ex: If (InRaster = 21, 22, 23, 24) and (inRaster2 = 11) - return 11) could anyone suggest the correct syntax for this kind of multiple conditions.


texts: Converting a large square polygon to smaller polygons

texts:
I have a large square polygon and I want to divide smaller polygons this polygon.
I have a large dataset that consists of points inside the polygon.
I want to calculate how many points inside the smaller polygons.
So how can I create a grid and export each gird polygons?
Here is the visualisation on the Google Earth.
Large green square is large polygon polygon. I want to create smaller square polygons like below blue squares and finally my goal is that to calculate how many red point in these blue polygons and sorted largest to smallest.
Below coordinates belong to green polygon.
<coordinates>
    29.31234019515055,40.93725819383494,0 29.34257575548109,41.03265424239162,0 29.20971921411221,41.05576636736616,0 29.17179939586544,40.96591648323777,0 29.31234019515055,40.93725819383494,0
</coordinates>


texts: Iterating and matching feature classes as an input by their names using ModelBuilder?

texts:
I am trying to iterate a series of feature classes (polygons) from two different/separated folders in model builder to convert a series of polygons using union and then the output for aggregate polygon tool. However, I am having an issue with the input file names. I have to match each entry with part of their names which are sharing a map#. In example both files are sharing map number like orange001 and apple001 or onion002 and banana002.
I need help on matching the name of each file entry and reading it as input for my model to create a desire output. I used 2 models as subtype and type which I can use the subtype and read the 1st file from my first folder and then I have to match its name with 2nd entry and grab the shp file and continue my work.


texts: Drawing specific layers of mapfile with MapScript?

texts:
Is it possible to call a draw on a mapFile with a list of layers to include with mapScript?
I have a bunch of layers defined on my mapfile but I don't want to call all of them. I don't use "STATUS DEFAULT", i would like to call specifics layers like we can do with the get parameters, for example: "LAYERS=LAYER1,LAYER2,LAYER3&...."

For now this is my code.

$oMapFinal = ms_newMapObj($filename) or die("die\n");
$oMapImageFinal = $oMapFinal->draw();


texts: Best way to serve vector tiles from data that changes frequently?

texts:
I'm working on a project where I process national radar data and generate a GeoJSON with the data. The data is generated every 5 minutes and I delete the oldest data before I add the latest data. I'm currently trying to decide the best way to go about serving vector tiles from this data. My current two solutions are:

Every 5 minutes I generate a new GeoJSON for the most current data, then using PostGIS I delete the oldest data from the table and then I insert the newest data. Once generated I serve tiles using a tool called Martin.
Every 5 minutes, I delete the oldest data from a GeoJSON stored in the file system, then I add the newest data to it. Then I use Tippecanoe to generate a .mbtiles file from the stored GeoJSON and serve the vector tiles using that.

From my research, I've found information for generating tiles on the fly and serving static tiles. But I don't need either. Right now, using PostGIS seems like the best option but it feels like overkill to generate the tiles every time they're requested. There isn't too much information out there for this scenario and I'm wondering if anyone could offer some insight into what might be the most optimal way to handle this.


texts: What causes ExecuteError: ERROR 000358: Invalid expression Failed to execute (SelectLayerByAttribute) on shapefile?

texts:
My script runs perfectly for feature classes, but not shapefiles. It crashes at the point of selection, provided there is a query.

gp.Workspace = "C:/someworkspace"
fc = "featureclass.shp"
query = "[SOME_TEXT_FIELD] = \"666\""
gp.MakeFeatureLayer(fc, "flyr")
gp.SelectLayerByAttribute("flyr", "NEW_SELECTION", query)

returns

  ExecuteError: ERROR 000358: Invalid expression Failed to execute
  (SelectLayerByAttribute).

Tested a similar simple example in the Python Shell. This runs fine with a geodatabase and feature class. I assume there is some sort of special treatment for shapefiles I'm not aware of, but arcgis desktop help does not specify anything.


texts: Is there a way to use File Geodatabases and Feature Classes, not just Shapefiles, in AutoCAD Civil 3D 2011 and any newer versions (i.e. 2013 and 2015)?


texts: Network Analysis - Service Area Tool doesn't consider all lines?

texts:
using the Service Area Tool from the Network Analysis (Toolbox) gives strange results. As you can see on the screenshot the calculated blue service area (network: streets; start points: points; path type: shortest; no advanced parameters set) does contain some of the side roads some others not although they should be clearly inside the service area. How is this possibile?

I already checked the streets layer (no gaps) and tried to run it with a multiline- and singleline-Layer (same result). Layers are in the same projected coordinate system. I'm using QGIS 3.6.3-Noosa.

Would be very happy if someone had an idea how to handle this. Thanks in advance!


texts:
I am trying to create a viewshed using the graphical modeler, the first part involves creating the viewpoints which works fine, however when I create the viewshed I get the following error. Also, when I do the exact same process outside of graphical modeler it works fine.

Traceback (most recent call last):

  File
  "C:/Users/Ben/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\viewshed_raster.py",
  line 245, in processAlgorithm interpolate = precision > 0) File
  "C:/Users/Ben/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\modules\visibility.py",
  line 178, in viewshed_raster data /= distance_matrix #all one line =
  (data -z - mxcurv) /mx_dist FloatingPointError: divide by zero
  encountered in true_divide

  Error encountered while running Viewshed


texts: Merge ("dissolve") OSM LineStrings from GeoJSON with ST_UNION

texts:
I would like to merge the continuous LineStrings from a GeoJSON export of an OpenStreetMap relation (get from Overpass Turbo):

  "geometry": {
    "type": "LineString",
    "coordinates": [

In another question I'm told to use:

ogr2ogr -f GeoJSON -explodecollections output.geojson input.geojson
-dialect sqlite -sql "SELECT LineString, ST_Union(geometry) as geometry FROM OGRGeoJSON GROUP BY LineString"

but since GDAL 2.2 "OGRGeoJSON" is no longer working.

How can I now dissolve such features?


texts: Adding legend title in plot function of R

texts:
I use plot function to plot a raster in R. The plot function adds a legend as default. However, the default legend added does not have a title. For instance I need to add the unit such as meters.

I cannot add a title to legend. How can I do it in plot function?

plot(raster, main = "text", sub = "This is elevation", lwd =5, col =my.palette, legend = T )


texts: Dynamically rendering polygon, point, and line features

texts: I have a Feature Collection consisting of features having id, color_code, etc as attributes. I have to render this feature with the color in the color_code attribute. Is this possible and please suggest some links?


texts: Dissolve line based on same id and same value keeping the order from start to end

texts:
I have the solution for doing this in R but would like to avoid using other tools.
Is it possible to dissolve lines with the value sections/features while keeping the start to end order and avoid dissolving the same values into one feature?
If you can imagine values within the same id: values= 1,1,1,1,4,4,4,1,1,1,5,5,5
I would like to get a 1,4,1,5 outcome. Would it be possible to preserve the start and endpoint of each feature?
If there is not a possibility to use dissolve can this kind of calculation be done in a Calculator? Basically creating a new group id.
Where new id would be for each value: 1,1,1,1,2,2,2,3,3,3,4,4,4
I will be able to use this new id for a dissolve.


texts: 3D Rayshader Plot looks compressed compared to static 2D plot

texts:
Sorry if this is an extremely obvious question, however I'm not sure why it's happening.

I'm trying to plot the provinces of Canada as 3D map, however whenever I render them, they end up looking mis-shapen. In particular, the North/South axis just looks compressed, or E/W is stretched. I tried a different coordinate system as well, however the issue seems to persist. I feel like I'm missing something elementary, but unfortunately I have no idea what. Would be grateful if someone could shine some light on this.

Some example photos include:

BC - Lat/Long

or, BC - UTM

The data for the plot I got from here and the code to plot them follows.

library(rayshader)
library(sp)
library(raster)
library(scales)
library(stringr)
library(maptools)
library(dplyr)
library(rgdal)
library(sf)
library(av)

province_name = "British Columbia"
abbrev = "bc"
path= sprintf("F:\\Datasets\\Geographical\\Topographical\\Canada\\%s\\", gsub(" ", "", province_name))

list.files(path=path_nun, full.names=TRUE) %>%
  (function(x) x[str_detect(x, "tif")]) ->
  elevations

dems = list()
for(i in 1:length(elevations)) {
  print(i)
  dems[[i]] = raster::raster(elevations[[i]])
}

all = do.call(raster::merge, dems)

# Get Province Shapefile
can = getData("GADM", country="CANADA", level=1)
coast = can[match(toupper(province_name), toupper(can$NAME_1)),]

# Crop the elevation to be the extent of the coastlines
elevation.sub = raster::crop(all, extent(coast))
elevation.sub.masked = raster::mask(elevation.sub, coast)

# 2D PLOTS
plot(elevation.sub.masked)

projection_string = "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
elevation.sub.masked.utm = raster::projectRaster(elevation.sub.masked, crs=projection_string)

# Rayshader requires matricies
elmat = raster_to_matrix(elevation.sub.masked)

# Reduce size
small_mat_small = reduce_matrix_size(elmat, 0.1)

# 3D PLOTS
zfactor = 10
small_mat_small %>%
  sphere_shade(texture = "bw") %>%
  add_shadow(ray_shade(small_mat_small, zscale=zfactor, multicore = TRUE), 0.3) %>%
  add_shadow(ambient_shade(small_mat_small, zscale = zfactor, multicore = TRUE), 0.2) %>%
  plot_3d(small_mat_small, zscale = zfactor, zoom = 1, theta=0, fov=0,  shadowdepth=-100, windowsize = 1000)


texts: QGIS processing script: Move points by specified distance and add segment

texts:
I try to move points (they come as WGS84, but I'm not sure if it works with unprojected CRS) by a specified distance and connect original and new point with a line segment with the below script, but get an error (EDIT: is fixed now..):

#ERROR: 'module' object has no attribute 'getobject'  --> fixed syntax error!

from PyQt4.QtCore import *
from qgis.core import *
from processing.core.VectorWriter import VectorWriter

#====================================
##[User scripts]=group
##input=vector
##output=output vector
##x_add=number 100
##y_add=number 100
#====================================

inputLayer = processing.getObject(input)   # formerly getobject()
features = processing.features(inputLayer) # formerly getFeatures()

provider = inputLayer.dataProvider()
writer = VectorWriter(output, None, provider.fields(), WKBLineString , inputLayer.crs())

for ft in features:
    line_start = ft.geometry().asPoint()
    line_end = QgsPoint(line_start.x()+x_add, line_start.y()+y_add)
    outFeat = QgsFeature()
    outFeat.setGeometry(QgsGeometry.fromPolyline([line_start, line_end]))
    writer.addFeature(outFeat)

del writer

Furthermore I'd like to be able to move points by specifying the distance in map units (m), I'm not sure if this is what's actually happening here.


texts: Calculate geometry not working in ArcGIS Pro?

texts:
I have a set of points whose x and y coordinates I would like to add to the attribute table. The points are stored as a feature class in a geodatabase and I have ensured that they have the proper co-ordinate system (World Mercator).

I have tried using both Calculate Geometry Attribute and Calculate Field to add the data to the columns named X and Y but both just yield <Null> even though the tool runs successfully. For Calculate field, I am using the following Syntax:

X=!SHAPE.CENTROID.X!
Y=!SHAPE.CENTROID.Y!

Where am I going wrong?

I started testing at ArcGIS Pro 2.2.0 but have updated to 2.2.3 and the problem persists.


texts: Selecting features of layer present just in one feature of another layer?

texts:
I'm a student approaching to QGIS 2.0.1-Dufour. My goal is to know how many features of a layer are present just in one feature of another layer. I tried both with the tools "spatial query" and "intersect" (see the picture) but I can't select just one feature.

How do I select the number of the red dots that are present in just one area defined by the black line?


texts: Getting warning when opening NetCDF file using gdal?

texts:
I am using windows 7(32-bit) and want to warp NetCDF (nc) files and save the output as GeoTIFF.

First I run gdalinfo first to get the info of the subdataset I need and added NETCDF in the beginning(e.g. gdalinfo NETCDF:"A2015069035500.L2_LAC_OC.nc") but was unfortunate because a warning appeared:

  Warning 1: Variable has 0 dimension(s) - not supported. gdalinfo
  failed - unable to open 'netCDF:A2015069035500.L2_LAC_OC.nc'.

Command line input:

  gdalinfo A2015069035500.L2_LAC_OC.nc

NetCDF is present when I tried gdalinfo --formats.

Does anyone have any idea about this?


texts: How to install pgRouting on PostgreSQL 9.2 located at /opt/PostgreSQL/9.2/ on Ubuntu Server 12.10?

texts:
I have installed PostgreSQL 9.2 on Ubuntu Server 12.10 using one-click installer from the EnterpriseDB. Then using the StackBuilder I installed PostGIS.

Now I also want to install the pgRouting extension to my database. So far I have found instructions on the following pages:

http://pgrouting.org/download.html
http://docs.pgrouting.org/2.0/en/doc/src/installation/index.html#ubuntu-debian
http://workshop.pgrouting.org/chapters/installation.html#pgrouting

sudo add-apt-repository ppa:georepublic/pgrouting-unstable
sudo apt-get update
sudo apt-get install postgresql-9.1-pgrouting

But the last command try to install postgresql-9.1 from the apt repositories:

# sudo apt-get install postgresql-9.1-pgrouting
Reading package lists... Done
Building dependency tree
Reading state information... Done
The following extra packages will be installed:
  libboost-thread1.49.0 libcgal9 libgeos-3.3.3 libgeos-c1 libproj0 libqt4-dbus libqt4-declarative libqt4-designer
  libqt4-help libqt4-network libqt4-opengl libqt4-script libqt4-scripttools libqt4-sql libqt4-sql-mysql
  libqt4-sql-sqlite libqt4-svg libqt4-test libqt4-xml libqt4-xmlpatterns libqtcore4 libqtgui4 postgis
  postgresql-9.1 postgresql-9.1-postgis postgresql-common proj-data qdbus
Suggested packages:
  proj-bin libqt4-declarative-folderlistmodel libqt4-declarative-gestures libqt4-declarative-particles
  libqt4-declarative-shaders qt4-qmlviewer libqt4-dev qt4-qtconfig oidentd ident-server locales-all osm2pgrouting
  pgrouting-workshop
The following NEW packages will be installed:
  libboost-thread1.49.0 libcgal9 libgeos-3.3.3 libgeos-c1 libproj0 libqt4-opengl postgis postgresql-9.1
  postgresql-9.1-pgrouting postgresql-9.1-postgis postgresql-common proj-data
The following packages will be upgraded:
  libqt4-dbus libqt4-declarative libqt4-designer libqt4-help libqt4-network libqt4-script libqt4-scripttools
  libqt4-sql libqt4-sql-mysql libqt4-sql-sqlite libqt4-svg libqt4-test libqt4-xml libqt4-xmlpatterns libqtcore4
  libqtgui4 qdbus

I want to install pgRouting on my existing PostgreSQL which is located at:

/opt/PostgreSQL/9.2/

I have also downloaded the source code. Following these instructions:

mkdir build
cd build
cmake -DWITH_DD=ON ..
make
sudo make install

When I run

# cmake -DWITH_DD=ON ..
-- PGROUTING_GIT_TAG: pgrouting-2.0.0
-- PGROUTING_GIT_BUILD: 0
-- PGROUTING_GIT_HASH: f26831f
-- PGROUTING_GIT_BRANCH: master
-- PGROUTING_VERSION_REVISION_NAME: f26831f master
-- PGROUTING_VERSION_REVISION_HASH: f26831f
-- UNIX=1
-- WIN32=
-- POSTGRESQL_EXECUTABLE is POSTGRESQL_EXECUTABLE-NOTFOUND
-- POSTGRESQL_PG_CONFIG is POSTGRESQL_PG_CONFIG-NOTFOUND
-- PostgreSQL not found.
CMake Error at CMakeLists.txt:122 (message):
   Please check your PostgreSQL installation.

-- Configuring incomplete, errors occurred!

I can't figure out a way to tell cmake to locate my existing PostgreSQL installation.

Any suggestions are very welcome.


texts: Complexify geography polygons in PostGIS

texts:
Is it possible to modify polygons (essentially, make them more complex by inserting additional points) so that their naïve point-linking rendering would be a closer match to their respective geodesic rendering?
Especially with such cryptic explanation, a picture is worth a thousand words...

Our use case is that we perform point-in-polygon outside PostgreSQL, by exporting the polygons, then performing PIP. The problem is that our PIP algorithm doesn't handle geodesics, therefore we'd like to be able to extract polygons from PostgreSQL which are closer to geodesic counterparts.
Note: no, for a myriad of reasons, we're not interested in doing PIP from inside PostGIS.


texts: getFeatureInfo is working only if line of polygon is clicked, does not work inside polygon

texts:
I have issue with getFeatureInfo function. I can get features of polygon only if it is clicked on line (border) of polygon. By clicking on line, i get features of few polygons, and i want to click inside of polygon and get feature info of that polygon, nothing else.
I'm using Openlayers 2.13 and GeoServer. Where can i fix this?


texts: How to enter survey coordinates into parcel fabric traverse grid?

texts:
I have a survey that only gives me the coordinate points for a parcel I need to draw into the parcel fabric. I am used to seeing bearings and distances - how do I take these coordinate points and turn them into bearings and distances without teaching myself a bunch of trig? Alternately, is there a tool in the parcel fabric wherein I can enter the points instead of the bearing and distance to create the parcel?
On a side note, does anyone have any recommendations for books/websites/resources that would help me learn about basics of surveying (as a user) to help me with math stuff like this in the future? Preferably at a "For Dummies" or ELI5 level...


texts: How to load all items from a template with pyQGIS in QGIS 2.6

texts:
Since QGIS 2.6, the LoadFromTemplate method seems not work any more (or works differently?)

This is my code to add items in my composer from a template :

self.composition = self.composerView.composition()
file1=QtCore.QFile('ComposerTemplate.qpt')
doc=QtXml.QDomDocument()
doc.setContent(file1, False)
self.composition.loadFromTemplate(doc, substitutionMap=None, addUndoCommands =False)

Up to 2.4 version, this code worked great.
But since the last QGIS update, nothing is imported from the template. I do not get either error messages

I've readen in QGIS changelogs, that the "Load from template" button is replaced by a "softer" "Add items from template" one. But it's in the GUI...
In the API doc, the LoadFromTemplate method is still present. And should still work?

Has someone an idea about what's happening?


texts:
I have two rows with monthly rain run off per watershed as in the picture below.

I'd like to obtain the average runoff between may and august for basin 218187 and 218265 (note that basin 218187 will be 0 but the attribute table contains more rows). To do that I've created a new field but I'm struggling to write the code in the field calculator.


texts: ModelBuilder Error 000354 Name contains invalid characters, using input selected from drop downs?

texts:
I've made a simple model in ModelBuilder, which is failing at the first step saying the name has invalid characters.

The failing tool is Copy Features:

The error is:

Executing (Copy Features): CopyFeatures "Geofabric Surface Network - V2.x\AHGFCatchment" in_memory\temp123 # 0 0 0
Start Time: Wed Mar 29 17:01:20 2017
ERROR 000354: The name contains invalid characters
ERROR 000354: The name contains invalid characters
ERROR 000354: The name contains invalid characters
Failed to execute (Copy Features).
Failed at Wed Mar 29 17:01:20 2017 (Elapsed Time: 0.22 seconds)

The output name in_memory\temp123 doesn't use any strange characters or bad words.

The input name was selected from the dropdown of layers available in the MXD:

The three 0 in the tool parameters are defaults and recommended by Arc.

I'm running this tool with some data selected from two layers - the purpose is to extract just the selected data to continue working on it.

What should I try next?


texts: ArcGIS API for JavaScript Web Optimizer: 'Queued' forever

texts:
Although ArcGIS web app optimizer behind the scenes uses the same dojo build system, it runs at least a twice slower. But the most unpleasant thing is it puts a building task in a queue with the “Your build has been queued and will start when resources become available. Feel free to close this browser tab; you'll receive an email when the build completes.” message. It seems that the queue never comes to an end. After the task gets status 'Queued', the next build task will have the same status, even the first one was canceled.
In such way it is impossible to create any build because of the queue.
I tried to create a build using ‘ArcGIS web app optimizer’ on Friday morning. The task was put to queue again. On Monday I have to cancel it because it was still in queue.
Then the ‘ArcGIS web app optimizer’ put in queue each build task I tried to create.  Further it was impossible to create any build by ‘ArcGIS web app optimizer’, even the simplest.

Please, help to figure out why the 'Queued' status appears and how to cancel it to create a new build.


texts: IGeometryBridge.GetPoints throws NotImplementedException. Am I doing something wrong?

texts:
I am working on some .NET code that processes lots of polygon geometries at the vertex level. In order to reduce COM interop overhead, I want to call IGeometryBridge2.GetPoints once instead of calling IPointCollection.Point repeatedly:

// using ESRI.ArcGIS.Geometry;
static IPoint[] GetVertices(this IPolygon singlepartPolygon)
{
    var pointCollection = (IPointCollection4)singlepartPolygon;
    var points = new IPoint[pointCollection.PointCount];
    var geometryBridge = (IGeometryBridge2)new GeometryEnvironment();
    geometryBridge.GetPoints(pointCollection, 0, ref points);
    return points;
}

geometryBridge.GetPoints(…) throws a NotImplementedException.

(I have tried initialising the points array with empty Point instances before calling GetPoints, but this didn't make any difference. I tried geometryBridge.QueryPoints(…), too, and this failed with a different exception.)

Am I doing something wrong, or is this more likely to be a bug?


texts: How to export KML from points plotted in javascript Google Maps API

texts:
I want to export a KML file from points plotted in Javascript using the Google Maps API. Is there a way to do this?

I would like the user to be able to log in and see Google Maps, then plot points on the map and finally export those points in a KML format.

Thanks.


texts: OSM2PO V5.1 - Routing from mid line not working with Ovl routers

texts:
I am using osm2po v5.1.0

I have combined 6 country pbf extracts as such:

java -jar osm2po-core-5.1.0-signed.jar prefix=sk cmd=tj slovakia-latest.osm.pbf
java -jar osm2po-core-5.1.0-signed.jar prefix=at cmd=tj austria-latest.osm.pbf
java -jar osm2po-core-5.1.0-signed.jar prefix=cz cmd=tj czech-republic-latest.osm.pbf
java -jar osm2po-core-5.1.0-signed.jar prefix=hu cmd=tj hungary-latest.osm.pbf
java -jar osm2po-core-5.1.0-signed.jar prefix=pl cmd=tj poland-latest.osm.pbf
java -jar osm2po-core-5.1.0-signed.jar prefix=ua cmd=tj ukraine-latest.osm.pbf
java -jar osm2po-core-5.1.0-signed.jar prefix=slovakia-plus sk at cz hu pl ua cmd=m
java -jar osm2po-core-5.1.0-signed.jar prefix=slovakia-plus cmd=sg

I have uncommented all these Ovl routes:

router.4.class=de.cm.osm2po.routing.OvlRouter
router.5.class=de.cm.osm2po.routing.OvlEdgeRouter
router.6.class=de.cm.osm2po.routing.OvlPoiRouter

I have used this URL to try and route starting away from a vertex (mid-line), using routerClassID=4:

  http://localhost:8888/Osm2poService?format=geojson&routerClassId=4&cmd=fr&source=48.92456,20.85100&target=48.92439,20.69805&findShortestPath=false&ignoreRestrictions=false&ignoreOneWays=false&heuristicFactor=0.0&maxCost=0.0&debugExpansion=false&key=Key&value=Value&tsp=

However, my results look exactly the same as using routerClassID=0 (DefaultRouter)

The end points of the routes are the same, lying on a vertex:

Am I missing a configuration? What am I doing wrong?


texts: Iterating through erase features while Input feature == output feature using ArcGIS ModelBuilder?

texts:
I have a polygon which represents a large area. Smaller polygons should be used to erase pieces step by step from of the large area polygon.
For this purpose, I use the Model Builder in ArcGIS.

My setting (please have a look at the picture I attached)

Iterator goes through a workspace where 'erase features' are stored and pass them to the 'erase tool'.
'erase too'l take the passed features as 'erase feature'.
'Input Feature' should be 'output feature' at the same time, because the process is used to delete small pieces of area from the 'input polygon'.

Problem: The tool delivers the message: "Output Feature Class is same as Input Features".

I tried to adapt this advice in Overwrite the input feature with the output feature with ArcPy, but storing the 'input feature' in in_memory causes the same message: "Output Feature Class is same as Input Features".

Is there a simple solution?


texts: OpenCV array is totally different from raster values of my TIFF in QGIS?

texts:
So I have a a GeoTIFF that I need to binarize. OpenCV is of course the choice program for the job, but I am having major issues.
Below, you will see a typical range for RGBnIR values in QGIS of my image:

First, the GeoTIFF is in UInt16, which I know is not appreciated by OpenCV. However, I can't seem to convert it/scale properly in QGIS to an 8bit unsigned integer. When I do, and open it up in OpenCV, there is simply no data in the bumpy array >0.
Confusingly, when I send the UInt16 data to OpenCV, I get these values, which make no sense (they should be up to a maximum of just over 1000 but the maximum is just 10):

I just want to create a simple binary image, but can't do that unless the number scales line up here...


texts: Open multiple files using fileopendlg() function

texts: I want to open multiple files using fileopendlg() function, when I code it selects only one file at a time. Is there any function or procedure to select multiple files at a time?


texts: Downloading GeoTIFF files from Soilgrids 2.0 using Linux yields files without spatial information

texts:
I've been trying to download global GeoTIFF files from this Soilgrids database in Linux using this tutorial. The process seems to work, but when I try to use the files they have no associated spatial information and are essentially just pictures. I'm an undergrad who is new to coding and does not have access to ArcGIS Desktop.
Here's the output when I run the gdalinfo function on one of the faulty GeoTIFF files:
Driver: GTiff/GeoTIFF

Files: ocs_0-5cm_mean_5000.tif

Size is 7962, 2902

Coordinate System is:

PROJCRS["Interrupted_Goode_Homolosine",

BASEGEOGCRS["GCS_WGS_1984_ellipse",

    DATUM["World Geodetic System 1984",

        ELLIPSOID["WGS 84",6378137,298.257223563,

            LENGTHUNIT["metre",1]],

        ID["EPSG",6326]],

    PRIMEM["Greenwich",0,

        ANGLEUNIT["Degree",0.0174532925199433]]],

CONVERSION["unnamed",

    METHOD["Interrupted Goode Homolosine"]],

CS[Cartesian,2],

    AXIS["(E)",east,

        ORDER[1],

        LENGTHUNIT["metre",1,

            ID["EPSG",9001]]],

    AXIS["(N)",north,

        ORDER[2],

        LENGTHUNIT["metre",1,

            ID["EPSG",9001]]]]

Data axis to CRS axis mapping: 1,2

Origin = (-19949750.000000000000000,8361000.000000000000000)

Pixel Size = (5000.000000000000000,-5000.000000000000000)

Metadata:

  AREA_OR_POINT=Area

  Code_version=v2.0.0
              Covariates=clm_wcl_srdyrsum,clm_wcl_s04rad,clm_wcl_p10tot,veg_mod_eviyravg,clm_wcl_p09tot,veg_mod_nppy15,    mor_env_demm,clm_wcl_s08rad,clm_wcl_srdyrstd,clm_wcl_s07rad,clm_wcl_bio19,clm_wcl_p05tot,clm_wcl_bio17,mo    r_mrg_twi,clm_mod_lstd03std,clm_mod_lstd04std,clm_wcl_bio16,clm_wcl_bio08,clm_mod_lstd09std,clm_mod_lstd1    1std,clm_wcl_bio18,mor_mrg_vdp,clm_wcl_p04tot,clm_mod_lstd05std,clm_mod_lstd12std,clm_mod_lstd10std,clm_m    od_lstd02std,clm_wcl_p12tot,clm_mod_lstd01std,veg_mod_evient,veg_mod_evirng,veg_mod_evimax,veg_mod_evievn    ,luc_gfc_trely10,mor_mrg_vbf

Litter_layers=FALSE

Model=Quantile Regression Forests

Model_type=ranger

Mtry=15

Number_trees=200

Outputs_version=RUN03

WoSIS_version=Data stream 7

Image Structure Metadata:

COMPRESSION=DEFLATE

INTERLEAVE=BAND

Corner Coordinates:

Upper Left  (-19949750.000, 8361000.000)

Lower Left  (-19949750.000,-6149000.000)

Upper Right (19860250.000, 8361000.000)

Lower Right (19860250.000,-6149000.000)

Center      (  -44750.000, 1106000.000) (  0d51'53.55"W,  9d56' 7.32"N)

Band 1 Block=256x256 Type=Int16, ColorInterp=Gray

NoData Value=-32768


texts: Generating random points near a location

texts:
In the following solution, why is there a sqrt(u)? Since u is just a random number in [0, 1) doesn't getting the sqrt just give you another smaller random number?

To generate points uniformly, randomly, and independently within a circle of radius r around a location (x0, y0), start by generating two independent uniform random values u and v in the interval [0, 1). (This is what almost every random number generator provides you.) Compute

d = r * sqrt(u)

theta = 2 * Pi * v

x = d * cos(theta)

y = d * sin(theta)

The desired random point is at location (x+x0, y+y0).


texts: How to draw two or multiple lines from one location to another without overlapping the lines?

texts: I want to draw two or many curved lines from fixed start point to another fixed end point on map. So, I started with geodesic but it draws overlapping lines. Then I tried to draw trajectory but not very useful as points could be any where in the world. How can I achieve this? I am using python and javascript. The below image is something I am aiming for (map is not added but imagine as if it is there :)). Also, is there any way to add directional arrow? Please help.


texts: Checking if service pack is installed using ArcPy?

texts:
Is there a way to check with arcpy whether or not a specific ArcGIS for Desktop service pack is installed on the local machine?

More specifically, my script can only work with ArcGIS 10.1 with SP1 or higher versions. In case of earlier versions, I want to display a message to the user and stop the script before it crashes.


texts:
I'm having some issues with Laspy and reading in .laz files.
Previously I have had no problem with this, but it has stopped working.

Initially I had a problem with the following error message: "laspy.util.LaspyException: Laszip was not found on the system"

I eventually got rid of this error using the top comment from here (although I already had the whole LAStools package so I don't understand why I needed to add it again).

However, this now throws up a new error with the following message(s):

  File "<ipython-input-3-025be98bbf71>", line 1, in <module>
    inFile2 = File('/path/to/data/file/myfile.laz', mode='r')

  File "/path/anaconda3/lib/python3.7/site-packages/laspy/file.py", line 64, in __init__
    self.open()

  File "/path/anaconda3/lib/python3.7/site-packages/laspy/file.py", line 75, in open
    self._reader = base.Reader(self.filename, mode=self._mode)

  File "/path/anaconda3/lib/python3.7/site-packages/laspy/base.py", line 273, in __init__
    self.setup_read_write(vlrs,evlrs, read_only=True)

  File "/path/anaconda3/lib/python3.7/site-packages/laspy/base.py", line 321, in setup_read_write
    self.data_provider.point_map()

  File "/path/anaconda3/lib/python3.7/site-packages/laspy/base.py", line 176, in point_map
    count=self.manager.header.point_records_count)

TypeError: a bytes-like object is required, not 'FakeMmap'

I can't seem to find any help on this.

Note that I am using Laspy v1.6.0.

Also note that I do not want to just use .las rather than .laz

Any help?


texts:
I am new to GeoTools. I was trying the samples in http://docs.geotools.org/latest/userguide/tutorial/raster/image.html.
I am not able to run the ImageLab sample code. I know it would be some lack in understanding or a simple mistake.

I am getting the following exception:

Exception in thread "main" java.lang.UnsupportedOperationException: Trying to get a reader from an unknown format.
    at org.geotools.coverage.grid.io.UnknownFormat.getReader(UnknownFormat.java:62)
    at org.geotools.demo.ImageLab.displayLayers(ImageLab.java:93)
    at org.geotools.demo.ImageLab.getLayersAndDisplay(ImageLab.java:77)
    at org.geotools.demo.ImageLab.main(ImageLab.java:51)

What actually needs to be done?  Also, is it possible to use normal jpeg image or any other way around?

In the GeoTools website the mentioned bluemarble.tif is not available and I have been using normal jpg image.

In another website while seraching for more info on this, I read that it uses a world file. What is a world file and how it is related to this and how do I use it?

Please help, I have been stuck with this for a long time.

Thanks in advance.


texts: GeoServer slow rendering of line features with border

texts:
I've published some shapefiles including line features using GeoServer and styled them with simple stroke and label. To apply a thin border around the roads, I've used another 'FeatureTypeStyle' according to the border section of SLD cookbook as follow:

<!-- border -->
<se:FeatureTypeStyle>
    <se:Rule>
      <se:LineSymbolizer>
        <se:Stroke>
          <se:SvgParameter name="stroke">#d9d5cd</se:SvgParameter>
          <se:SvgParameter name="stroke-width">7</se:SvgParameter>
          <se:SvgParameter name="stroke-linejoin">round</se:SvgParameter>
          <se:SvgParameter name="stroke-linecap">round</se:SvgParameter>
        </se:Stroke>
      </se:LineSymbolizer>
    </se:Rule>
</se:FeatureTypeStyle>
<!-- line -->
<se:FeatureTypeStyle>
    <se:Rule>
      <se:LineSymbolizer>
        <se:Stroke>
          <se:SvgParameter name="stroke">#F5F1E6</se:SvgParameter>
          <se:SvgParameter name="stroke-width">5</se:SvgParameter>
          <se:SvgParameter name="stroke-linejoin">round</se:SvgParameter>
          <se:SvgParameter name="stroke-linecap">round</se:SvgParameter>
        </se:Stroke>
      </se:LineSymbolizer>
      <se:TextSymbolizer>
        <se:Label>
          <ogc:PropertyName>name</ogc:PropertyName>
        </se:Label>
        <se:LabelPlacement>
          <se:LinePlacement>
            <se:PerpendicularOffset> 0 </se:PerpendicularOffset>
          </se:LinePlacement>
        </se:LabelPlacement>
        <se:Fill>
          <se:SvgParameter name="fill">#523735</se:SvgParameter>
        </se:Fill>
        <se:VendorOption name="followLine">true</se:VendorOption>
        <se:VendorOption name="maxAngleDelta">90</se:VendorOption>
      </se:TextSymbolizer>
    </se:Rule>
</se:FeatureTypeStyle>

When I applied line border, the rendering time of WMS tiles increase in comparison with simple lines situation. How can I achive an optimal rendering performance?


texts: Create mapping between overlapping linestrings

texts:
What I have is two shapefiles, describing roads in a city. The shapefiles are made of linestrings, each covering a section of road.

The road segments in shapefile 1 are generally several times larger than the segments in shapefile 2.
The road segments break at road junctions for both shapefiles.
The shapefiles do not describe every road in the city.
Both shapefiles describe the same city, although there is some difference in the roads they cover. Some roads decribed by segments in shapefile 1 are not present in shapefile 2, and vice versa.

What I want is to create a many to one map between the shapefiles. I would like a table that lists each segment in shapefile 2, and it's closest match in shapefile 1. By closest match I mean closest in similarity, not necessarily in distance. In situations where a segment in shapefile 2 overlaps two segments in shapefile 1, the segment it overlaps the most would be chosen. Where there is no matching segment in shapefile 1, the table entry would be null.

I'm not sure how to best approach this task. My initial plan is to create a new polygon layer using the buffer tool in QGIS, from the features in shapefile 1, and then somehow find which buffer the linestrings from shapefile 2 predominantly fall in. I don't know if there are any tools which would help with this second step, or tackles the whole problem. I know how to use QGIS, and Python.


texts: how to calculate polygon area in arcengine

texts:
i'm new to use Arcengine 9.3 . i would like to create an application to calculate the selected polygon in combobox such that the user click to a button and polygon area saved in dbf .
how can i write the codes with c# in arcengine ?
thank you


texts: Calculation of focal statistics at the edges of the region

texts: When I use focal statistics setting the neighborhood of each processing cell to be a circle with 800m radius (map units, not number of cells), how are the statistics calculated in the edges of my region, when in that radius there are ΝoData pixels? In my case, I have selected for each neighborhood to be defined by a circle with 800m radius around each processing pixel. What I get in the edges is something like a "buffer" around my region that equals to the circle radius. Is this normal? How is it explained?


texts:
Can you do Raster Color Slices (like in ENVI) in QGIS?

I.e. Where the histogram can be divided up into discrete units to represent different classes based on the values of the pixels.


texts: How do I extract beginning/ending points for polyline shapes?

texts:
I need to extract the beginning/ending coordinates of every shape within a polyline shapefile. It would be best to actually add them as values within the shapefiles .dbf (I can live without doing this).

What is the easiest way to do this in ArcMap? I know that XTools has a plugin that does this really easy, but I don't want to purchase software for something that seems like it should be fairly trivial.

I am currently using ArcMap 10.


texts: Merge layers without losing excess observations in ArcGIS?

texts:
I am currently using Arcgis to map UK Election Data. I merged a shapefile and a datatset together and have starting looking at trends. See image below.

Then I looked at the attribute table of the merged item and realised I'd done something very stupid; My dataset had 20,000 observations, whereas my district map of the UK contained just over 400 items, since it's based on local authority level. In looking at the attribute table I realised that 400 records have been merged and the rest of the observations seemingly not there.

Thinking about it now, of course this would happen. I suppose what I wanted instead is a density or average of all the records relating to a particular polygon. For example, every observation contains a 1 - 10 score on a question around how succesfull the government has been. Obviously I don't want tp map just one observation per area, but instead want to aggregate ALL OBSERVATIONS per area.

Therefore my question is how to accomplish this in ArcGis, and how to merge layers without deleting all the 'excess' observations unneeded by the shapefile?


texts: ArcCatalog not showing view just created?

texts:
I have a geodatabase with a feature layer of parcels and a scalar database of consumption data for those parcels. Using SQL Server 2012r2 and ArcGIS 10.5.

In ArcCatalog, I open the parcel database connection and create a view which queries the consumption data.  The view creation completes without error, but ArcCatalog does not show the view (as a "table") in the connection.  Refreshing the connection, disconnecting and reconnecting the DB connection doesn't help.

Yet, if I try to create the view again, I get error "table already exists", and if I open the geodatabase from SQL Server manager, I can see the view and even run it and see expected data.

Why won't ArcCatalog show the view?

I suspect it doesn't enumerate DB schema from the geodatabase that it doesn't have SQL level access to, but connection to the geodatabase has full owner privileges (using SQL rather than Windows authentication).

Once the view is visible in the geodatabase, I plan to use it to define a relationship class with which I hope to be able to visualize consumption by (groups of) parcels, showing numbers and bars on a map, or color-coding parcels by consumption.


texts: WorkspaceFactory not defined in Python comtypes for ArcObjects?

texts:
I'm trying to programmatically export a symbolized raster to an RGB TIF using ArcGIS, i.e.: perform the equivalent of Export Raster Data in ArcMap with Use Renderer and Force RGB activated.

I followed the beginning of Use python to disable "Scale map elements proportionally to changes in page size" to access ArcObjects in Python.

I'm trying to re-create the code in this discussion to eventually create a RasterLayerExport object. Before that, though, I find myself incapable of creating a WorkspaceFactory.

I tried running:

import arcpy

from snippets102 import *
from comtypes.client import GetModule, CreateObject

import comtypes.gen.esriFramework as esriFramework
import comtypes.gen.esriArcMapUI as esriArcMapUI
import comtypes.gen.esriCarto as esriCarto
import comtypes.gen.esriDisplay as esriDisplay
import comtypes.gen.esriGeoDatabase as esriGeoDatabase

wkspf = CreateObject(esriGeoDatabase.WorkspaceFactory, interface = esriGeoDatabase.IWorkspaceFactory)

But this gives me the following error:

Traceback (most recent call last):
  File "C:/Users/AndyPL/.PyCharmCE2017.2/config/scratches/scratch_1.py", line 14, in <module>
    wkspf = CreateObject(esriGeoDatabase.WorkspaceFactory, interface=esriGeoDatabase.IWorkspaceFactory)
  File "C:\Python27\ArcGIS10.4\lib\site-packages\comtypes\client\__init__.py", line 238, in CreateObject
    obj = comtypes.CoCreateInstance(clsid, clsctx=clsctx, interface=interface)
  File "C:\Python27\ArcGIS10.4\lib\site-packages\comtypes\__init__.py", line 1225, in CoCreateInstance
    _ole32.CoCreateInstance(byref(clsid), punkouter, clsctx, byref(iid), byref(p))
  File "_ctypes/callproc.c", line 945, in GetResult
WindowsError: [Error -2147221164] Class not registered

Odd, since creating a PlugInWorkspaceFactory works fine. Does anyone else get this error?


texts: Counting values of each classification code using QGIS

texts:
I feel like this is an easy fix, but I'm struggling to just get a count of how many attributes I have for each value in a certain column for an Attribute table.
I have a column of landfire classification code is associated with my GPS points, and I want to know how many points I have that are classified as 7610, how many points are classified as 7611, and so on for each unique code.
Anyone know or have ideas?
I'm using QGIS 3.16.


texts:
I have layer with a lot of points being displayed.
It does not make sense to display these points until the zoom level is quite close.

How can I configure this layer to not show any points until a certain zoom level has been reached?

I am using geoserver 2.10 with the CSS extension installed.


texts:
ee.Image([1,2])

evaluates to two band constants image.

var myPoint = ee.Geometry.Point([13.683294059500382, 51.438313436248876]);
var coords = myPoint.coordinates()

evaluates to list of coordinates JSON
but

var constImage = ee.Image(coords)

print(constImage)

evaluates to list and not to two band image. Why? How can I create a constant multiband image from a ee.Geometry.coordinates() call?

EDIT: the background is that I want to calculate an image of euclidian direction from a defined point. My solution is to get x y images using ee.Image.pixelCoordinates and apply something like Math.atan(x-cx, y-cy), where cx and cy are the coordinate of myPoint. Thus, I need a constant images of cx and cy. The rest is then just band math and dealing with projections. I could just copy paste the coordinates, e.g. ee.Image(cx).addBands(ee.Image(cy)), but the goal is to create an app in which the center point is defined by a click. So I want to be able to create an ee.Image from myPoint.coordinates.

var geometry = ee.Geometry.Point([13.683294059500382, 51.438313436248876])
var proj = 'EPSG:32633'
var geometry = geometry.transform(proj)

//var constImage = ee.Image(geometry.coordinates())

var centerImage = ee.Image(408482.29).addBands(ee.Image(5699391.17)).rename('cx', 'cy')

var latlonImage = ee.Image.pixelCoordinates(proj)
var diffImage = latlonImage.subtract(centerImage)
var theta = diffImage.select('y').atan2(diffImage.select('x')).multiply(180).divide(3.14)
Map.addLayer(theta)


texts: PostGIS Clustering. Group by geographic point location proximity

texts:
I have a table where i keep user's location history. The system gets user's current location each minute and writes it into "user_locations" table.

Later on; when showing the locations that the user has been; this is what it seems like when selecting all records from "user_locations":

Apparently the user has been too lazy all day but it ends up so that i have to handle the outcomes of his laziness -which he is unaware of i think-.

Okay. What i want to do is to show only one marker in a (let's say) 50 meter radius and the marker label to indicate: "The user has been here from 02:00 to 23:00". Of course the dates would be the min and max date values of the records in the cluster.

How can i achieve this?

Here is my table structure:

Table Name: user_locations
Columns: id, user_id, created_at (timestamp), location (Point, geography, srid:4326)

One important thing to note is to break the cluster once the user has gone out of the radius. I mean; let's say the user has been sitting in his office from 08:00 to 13:00 and went to lunch at 13:00 and came back to office at 14:00 and took out of the office at 18:00. This would mean 3 clusters. 2 on office and 1 on lunch restaurant. So the clustering would be based on location but a location outside of the given perimeter would break the clustering and even if later the location becomes in the same place, it would start another cluster.

I hope i was able to tell what is on my mind.

Regards.


texts: How to optimize a WMS request by limiting number of tiles?

texts:
I want to display the layer using WMS and meanwhile I want to limit the number of tiles with their respecive BBox when map is loaded. At the backend I am using geoserver.

I can limit the number of features to load for WFS using "Maximum number of features" option in WFS setting in geoserver. But I am somehow lost to find the way of limiting number of tiles in WMS.
If I do like:

var bound = new OpenLayers.Bounds(6.9228, 50.9251, 6.9911, 50.9574);

and if i go and play around with different lat-lon bbox i can see in tomcat log file that, get request for tiles with respective BBOX are varying depending upon how deep i am looking for.

Is there any config setting in geoserver or openlayers api to control the tiles?
Thanks


texts: How can I limit a Hive spatial binning query to a specific map extent?

texts:
I'm trying to do a HQL query based off a map extent. I cant get this to work, what am I doing wrong? I'm using the taxi_demo in the ESRI example.

SELECT ST_Bin(0.001, ST_Point(dropoff_longitude,dropoff_latitude)) bin_id, *FROM taxi_demo WHERE ST_Contains(ST_Envelope(-8238484.195726455, 4969891.9664984355, -8235431.491909235, 4971962.932233592), ST_Point(dropoff_longitude,dropoff_latitude))

Am I doing something wrong or fat fingering something? It says: No matching method for class com.esri.hadoop.hive.ST_Envelope with (double, double, double, double). Possible choices: FUNC(binary)


texts:
I have a project with a bunch of large, high-resolution Tiff images that I will tilling with gdal2tiles and use on Openlayers.
Those Tiff images have been taken on a drone. Example:

While I tilling this by gdal2tiles and use Openlayers to show up, the resolution is low and the quality is very bad with the same zoom level, different from my expectations.

The image layers look like blur, something makes visible mask and the detail is gone. I use gdal on OS Shell and the code is:
gdal2tiles -p raster --tilesize=256 abc.tif output

How can I improve the quality and detail of this case?


texts:
I have gotton some help in creating a custom control for leaflet, but I am still having trouble upgrading to the CSS styles used by the latest versions of Leaflet. Basically, all I need is a icon in the white box, with the drop shadow around it.

Here is what I have working so far, take a look and feel free to fork it:

http://codepen.io/DrYSG/pen/zJsiG


texts: Any idea if it is possible to have an esri-asci file with a cellSize with 2 values specified? From what I read is not, but I have one file that is like that and I wader how should I import it in my postGis raster table.


texts: Find the point distance along a streams polyline?

texts:
Working in ArcGIS.

I am looking for a way to find the distance of a point along a streams polyline. I read Finding distance of point along polyline? and it seems to be similar to what I need...the distance of something downstream.

I went ahead and clipped my stream lines to only downstream to avoid any confusion on my part. I also created start points of each stream polyline.

I tried to use the 'Locate Features Along Route';  I get the distance the point is from the line...not how far down it is from a starting point...

I do not have the Network Analyst License.

Example of dissolved streamlines...gives not the closest route

Example of odd distances after minimizing available routes. Is point 1 really 81 feet down...if point 2 is only 50?


texts: Naming output same as input in ModelBuilder?

texts:
I have a model which has 6 select queries for different feet of contours.

I want to be able to give the same name as the input town/city for the output.

How can I do this?


texts: QgsVectorDataProvider capabilities not available in QGIS

texts:
I am trying to use deleteAttributes on my vector layer to delete attributes, but it does not work the attributes are not deleted.

layer.dataProvider().capabilitiesString() gives the following output:

u'Create Spatial Index, Fast Access to Features at ID, Simplify Geometries,
Simplify Geometries with topological validation'

which confirm deleteAttribute is not available. What are the reasons why the function (and the others) are not supported? I am using QgsVectorDataProvider.

The layer is created by joining an attribute table (table) (imported from a CSV) to a vector layer ("layer1") using

res= processing.runalg("qgis:joinattributestable", layer1, table, "ID", "ID", None)
layer = QgsVectorLayer(res['OUTPUT_LAYER'],"join_layer","ogr")

"layer1" has all capabilities in QgsVectordataProvider

u'Add Features, Delete Features, Change Attribute Values, Add Attributes,
Delete Attributes, Create Spatial Index, Fast Access to Features at ID,
Change Geometries, Simplify Geometries, Simplify Geometries with topological validation'


texts: Merging raw data from Landsat 5, Landsat 7 and Landsat 8

texts:
This question concerns the merging of the following Google Earth Engine collections:

USGS Landsat 5 TM Collection 1 Tier 1 Raw Scenes (LANDSAT/LT05/C01/T1)
USGS Landsat 7 Collection 1 Tier 1 Raw Scenes (LANDSAT/LE07/C01/T1)
USGS Landsat 8 Collection 1 Tier 1 Raw Scenes (LANDSAT/LC08/C01/T1)

I want to extend the "Landsat Simple Composite" example by loading similar collections from Landsat 5 and Landsat7, so that more years' composites are available.  I tried to first merge these three collections into a big one, before passing it to the Simple Composite algorithm:

// Composite 6 months of Landsat 8.

var l5 = ee.ImageCollection("LANDSAT/LT05/C01/T1");

var l7 = ee.ImageCollection("LANDSAT/LE07/C01/T1");

// Note that the input to simpleComposite is raw data.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1');

var lall = ee.ImageCollection(l5.merge(l7).merge(l8));

// The asFloat parameter gives floating-point TOA output instead of
// the UINT8 outputs of the default simpleComposite().
var composite = ee.Algorithms.Landsat.simpleComposite({
  collection: lall.filterDate('2015-1-1', '2015-7-1'),
  asFloat: true
});

// Pick a spot with lots of clouds.
Map.setCenter(-47.6735, -0.6344, 12);
// Display a composite with a band combination chosen from:
// https://landsat.usgs.gov/how-do-landsat-8-band-combinations-differ-landsat-7-or-landsat-5-satellite-data
Map.addLayer(composite, {bands: ['B6', 'B5', 'B4'], max: [0.3, 0.4, 0.3]});

However, this returned the error:

Layer 1: Layer error: Image.visualize: No band named 'B6'. Available band names: [B1, B2, B3, B4, B5, B6_VCID_1, B6_VCID_2, B7, B8].

I see that the bands from these three collections are not the same, but on the other hand, l5.merge(l8) runs with no errors. What's the best way to merge these collections?


texts: Creating a snake-like labelling pattern with sequential numbers (ArcPy)

texts:
I am new to Python and ArcPy.
I've got rows of points that go up in sequential number starting from the left to right. I need every second row to be flipped so that it goes right to left. Every point has an ID field denoting what row it belongs to (ORIG_FID). I can't figure out how to select that value in order to flip it. I'll also need to get the highest value in that particular row and plug it into pStart. Here is what I have:

def unique_values(table, field):
with arcpy.da.SearchCursor(table, [field]) as cursor:
    return sorted([{row[0] for row in cursor}])

values = unique_values(cruise_points, "ORIG_FID")

with arcpy.da.SearchCursor(cruise_points, "ORIG_FID") as cursor:
    for row in cursor:
        if (row[] % 2) != 0:
            select = arcpy.SelectLayerByAttribute_management("cruiselyr", "NEW_SELECTION", ' "ORIG_FID" = unique[]')
            codeblock = ('\n'
                         'rec = 0\n'
                         'def autoIncrement():\n'
                         '    global rec\n'
                         '    pStart = 12\n'  # HIGHEST VALUE IN THE ROW
                         '    pInterval = -1\n'
                         '    if (rec == 0):\n'
                         '        rec = pStart\n'
                         '    else:\n'
                         '        rec += pInterval\n'
                         '    return rec ')
            arcpy.CalculateField_management(select, "Label", "autoIncrement()", "PYTHON_9.3", codeblock)


texts: Is there an upgrade for Midvatten that can handle Mac OSX High Sierra. Sierra had no issues loading, but new installation of QGIS 2.18.15-4 bounces the plug-in.


texts: Point Feature Percentiles: National Rank vs Local Rank in ArcGIS

texts:
I have an idea in my head, but need assistance in how to process my data to get the answers I need.

Restrictions: Where I work we use Esri products, and third-party packages and add-ins are prohibited, so packages such as R or QGIS are not available to me. Python or ModelBuilder are available.

Problem Statement: I have a nation-wide point data set with values (approximately 60k points). The higher the value, the more "important" the site. Nationally, a site may rank low, but I am interested in how that site ranks to those nearby it. The goal is to find the local "big dog" sites.

Minimum value is 0. In practice, the maximum value is 119, but the theoretical maximum value is unlimited.

I do not want to divide the analysis by state, county, or any other political boundary. That would provide a scenario whereby a site would ignore a nearby site just over a border. The question I am trying to answer is "How do I compare to all of the sites within X miles of me?"

I would like to have three columns in the table:

1) Value:
This is the existing value

2) National Percentile:
A normalized value between 0-100 that indicates rank nationwide
ex: 99; 50; 23

3) Local Percentile:
A normalized value between 0-100 that indicates rank locally
Ideally, I would like an assessment done for each point that selects all points within X miles and determines the percentile of that point for the set of points within X miles of that one. Let's assume X = 100 miles
ex: 99; 50; 23


texts: Using GDAL to merge 3 bands to form a multispectral image?

texts:
I have 3 separate *.tif image files as grayscale. They are representing data acquired with a scanner. Each of them is a separate band of the scanner; red, green, and blue. I would like to reassemble them in a single RGB *.tif image file.

Using gdal_merge.py -o RGB.tif -co PHOTOMETRIC=RGB red.tif green.tif blue.tif

doesn't work; I end up with a single grayscale image which size is the sum of all individual images.
How to achieve this simple operation?

(They are not geotiff, just simple tif files; there is no tfw)


texts: Resolve getting a sorted list when using 'getFeatures()' method

texts:
I have this code which allows to insert features in a list:

liste = []
for b in layer2.getFeatures():
    for a in layer1.getFeatures():
        if a.geometry().intersects(b.geometry()):
            liste.append(a[12])

print liste

Normaly I expect this result [67,34,35,68],
but the output of my code give me a sorted list [34,35,67,68].
How can I resolve this problem?

THE COMPLETE CODE:

sel=[]

for a in layer1.getFeatures():
    for b in layer2.getFeatures():
        if a.geometry().intersects(b.geometry()):
            sel.append(a)
            print a['CODE_CH1']

b = len(sel)
d = QgsDistanceArea()
d.setEllipsoid('WGS84')
d.setEllipsoidalMode(True)

for i in range(b-1):
    j=i+1
    if i<j:
        a = d.measureLine(sel[i].geometry().asPoint(),sel[j].geometry().asPoint())
        print sel[i][12],sel[j][12],a


texts:
I have this data that is a result of intersections and areal calculations. In the attribute table drawn below, you can see the column Kelurahan has many duplicates because the mean prices differ based on the area. So I calculated the price_weights to measure price against proportion of area.

How do I simplify this, like add all the price_weights for each specific Kelurahan and then have only single Kelurahans?

I want to have one string price_weight for Kelurahans. So there is one row of Ubo-ubo for instance and the price_weight value is 61048,775+141766,763+......

Please leave a comment if something is unclear.


texts: Displaying message error if field is null using QGIS

texts:
I'm creating a form on some features in GGIS.

How do I display a message error if a field of a feature is null?

I do not want to use QtCreator.


texts: QGIS' 3.x - Atlas feature not printing second pdf page

texts:
So, I have created a layout managed by Atlas; every layout has 2 A4 pages, one with the atlas view and another with a table based on visible objects. It works perfect when printing to paper (duplex), but when I try to export as PDF I only get the first page, so the tables are lost.
I want to print a single PDF with both pages, for every atlas object
(also, I dont mind exporting atlas objects one by one, but I need them to have both pages)


texts: Converting shape Polyline-M to Polyline using ArcGIS Desktop?

texts:
I have some pretty old software that I am forced to use (no option). It's old enough that it doesn't recognize shapefiles of type Polyline-M, it only recognizes shapefiles of Polyline.

I have a shapefile that is completely Polyline-M that I need to convert to plain old Polyline.

Does anyone know how to do this in ArcMap 10?


texts: I'm just wondering if there's any way to easily lock the map scale when panning etc. For example, suppose I set the scale to 1:100,000 - panning with the mouse will then alter that slightly, which can be a real nuisance sometimes.


texts: Applying curve using ArcGIS API for JavaScript?

texts:
In desktop, there is http://pro.arcgis.com/en/pro-app/help/editing/create-arc-segments-and-curves.htm to create arc and curves, but I can't find any tutorial about how to make and create arc or curves on ArcGIS API for Javascript.

Is there any way to do that?


texts: Displaying 3rd party FeatureSet geometries on jupyter notebook webmap using ArcGIS Python API

texts:
How can I add all geometries included in a FeatureSet to a web map?

Trying to add geometries (i.e., shapes) to a Jupyter notebook web map (see code sample below).

Adding features from the ArcGIS server is a breeze. But what about adding from third parties?

My geometries are first imported as a .geojson, which I convert to an arcgis.feature.FeatureSet.

For each feature in the set, I'd like to use the coordinates of each shape to add that shape to a web map.

Setup: latest anaconda python with arcgis API v.1.4 on a Linux 16.04 box.

import arcgis
from IPython.display import display
import requests

# link build
HI_geo_stem = 'https://opendata.arcgis.com/datasets/'
HI_parcels_json = '1eb5fa03038d49cba930096ea67194e0_5'
parcels_link = ''.join([HI_geo_stem, HI_parcels_json, '.geojson'])

# read in parcel json
raw_parcels_json = requests.get(parcels_link).json()

# convert json to arcgis FeatureSet
hi_arc_parcels = arcgis.features.FeatureSet.from_geojson(raw_parcels_json)

# display a web mapgis = arcgis.GIS("https://www.arcgis.com")
m = gis.map("Hawaii", zoomlevel = 10)

# now display hi_arc_parcels in web map m ??


texts: Looking GeoServer layer name specification

texts:
I am using GeoServer 2.18.0 along with PostgreSQL database.
My most of the data is coming from the database. Certain table name starts with special character like "%", few table names contains "&" in between as layer name.
Services using database tables : WMS, WFS & CSW
Is there any layer naming convention document available for GeoServer? So that I will follow best practices.
Found following :
http://osgeo-org.1560.x6.nabble.com/Re-Geoserver-devel-wms-getCapabilities-special-characters-td5160099.html
https://docs.geoserver.org/latest/en/api/#1.0.0/layers.yaml
Styling WMS layers in GeoServer with fields using "%" character


texts:
I have a feature class that has been processed and says that the projection is in WGS 1984 Web Mercator Auxiliary Sphere.

However, I suspect that it is NOT truly in WGS 1984 Web Mercator Auxiliary Sphere. I suspect that it is still in its pre-processed projection, GCS North American 1927. Unfortunately I cannot get the source data prior to being processed and still in its original format of GCS North American 1927.

How can I accurately test to see if this dataset is truly in WGS 1984 Web Mercator Auxiliary Sphere or GCS North American 1927?

ArcGIS Desktop 10.3.1


texts:
Similar to the question on how to extract highways from osm data I posted before, I now have a new one: I want to extract all features tagged with railway=rail from an osm.pbf file using ogr2ogr. In contrary to the highway key, which is a header in the attribute table, the railway is not. I tried to change this by editing the osmconf.ini file, but I always get the error "ERROR 1: "railway" not recognised as an available field."
ogr2ogr -where "railway='rail'" CA-rails.shp central-america-latest.osm.pbf lines

Even though, the relevant line in the [lines] part of osmconf.ini was changed.
# keys to report as OGR fields
attributes=railway

So, the question arises here for me: Did I still not find the correct osmconf.ini file, and if so, where is the related file for OSGeo4W saved. Currently, I only find osmconf.ini files installed for R.
Otherwise, how can I extract the railway key which is saved under "other_tags" using the command line?
The line I used so far was:
ogr2ogr -where "other_tags='railway=rail'" CA-rails.shp central-america-latest.osm.pbf lines


texts: How to set flags when re-use QgsFeatureRequest object

texts:
I have a standalone PyQGis application that is working fine. When I query a vector layer via QgsFeatureRequest class, I would like to instantiate one object and reuse it. But I am not sure how to set- and unset- some flags. For example,

request = QgsFeatureRequest()
expr = QgsExpression(expr_text)
request.setFilterExpression(expr_text)
if no_geometry:
   request.setFlags(QgsFeatureRequest.NoGeometry)
else:
   # does the following line undo the above 'NoGeometry' setting?
   request.setFlags(QgsFeatureRequest.NoFlags)

if subset_attribute:
   request.setSubsetOfAttributes([ind])
else:
   # here, what should I do to undo the above 'setSubsetOfAttributes'???

I guess the bigger question is whether I should keep creating new implicit QgsFeatureRequest objects when query or should I keep one instance and change its flags as above, which way is better?


texts: Projecting many raster files in folder and then putting them into ArcGIS Pro project using ArcPy?

texts:
I want to project many raster files in a folder (but they appear in the content of my ArcGIS project). Here is the screenshot.

These raster files have an unknown coordinate system. When I try to project them manually, the input and output coordinate system is shown below. And I have successfully project several raster files manually.

As you can see, the input and output coordinate system is the same "GCS_WGS_1984". I can manually project these raster files but want to find how to accomplish this using python. Also, after projecting, how to make sure they will appear in my gdb project for future use.

I tried to use the following code but it did not work.

import arcpy
arcpy.env.workspace = r"D:\STAR\info_data_test"
#this is the folder where raster files are
arcpy.env.overwriteOutput = True #overwrite the exsiting files

rasterfiles = list(set(arcpy.ListDatasets("b*","Raster")))
#list all the raster files inside the folder with name starting "b"

base = "D:/STAR/assignment/assignment7/Zhenyu_Assign7_GISpro/Zhenyu_Assign7_GISpro.gdb/b90w90r65s400Pro"
#using the exsiting projected raster file as base
out_coor_system = arcpy.Describe(base).spatialReference
#get the coordinate system of this base raster file

for raster in rasterfiles:
    outraster = raster + "Pro"
    arcpy.ProjectRaster_management(raster,outraster,base)

I do not know how to project these raster files and then put them into my gdb project for future use.


texts: Different point buffer length left and right in Arcmap

texts: How am I supposed to split my buffer at different length from left and right. I need to make buffer 100m left and 300 m on right. Green is 100m while blue is 300m.


texts:
I am running GeoServer 2.4.3 on Ubuntu Server 12.04 and tomcat7 and I want to access it on mydomain.com. I have searched a lot and found the all the solutions suggest to use a ProxyPass. One of the posts on GIS Stack Exchange that has a similar problem is this one but the website to the solution provided in the post is down!

This is what I have been doing so far:

Added the follwing to /etc/apache2/sites-available/default

 ServerName mydomain.com
 ProxyRequests Off
 ProxyPreserveHost On

 <Proxy *>
     Order deny,allow
     Allow from all
 </Proxy>
 ProxyPass /geoserver http://servername:8080/geoserver
 ProxyPassReverse /geoserver http://servername:8080/geoserver

Then I have edited /var/lib/tomcat7/webapps/geoserver/WEB-INF/web.xml

<context-param>
      <param-name>PROXY_BASE_URL</param-name>
      <param-value>http://mydomain.com:8080/geoserver/</param-value>
    </context-param>

Finally, I have added In GeoServer > Global Settings > Proxy Base URL >

http://mydomain.com:8080/geoserver/

Now, when I am trying to enter mydomain.com:8080/geoserver I get Unable to connect error!

What I have been doing wrong? Any ideas?


texts: Raster brightness, contrast, saturation in PyQGIS

texts:
I'm trying to find a way to change these properties (especially contrast and brightnes) and I've found a QgsBrightnessContrastFilter class, but I can't fugure out how to apply it to my raster layer.

I think that can be done with:

contrastFilter = QgsBrightnessContrastFilter()
contrastFilter.setContrast(100)

But what now? How to set this filter to my layer?