texts:
I am running the Extract Values to Table tool from the Geostatistical Analysis Toolbox in ArcGIS Desktop 10.6.1 on Windows 10.
My inputs are a polygon layer with ~70.000 polygons and a raster layer with around 612,160,000 cells.

  Executing: ExtractValuesToTable HULL_polygon DGM10_slope
  J:\XXX\HULL_slope.dbf
  DO_NOT_ADD_WARNING_FIELD
    Start Time: Mon Sep 30 08:18:12 2019
    ERROR 999998: Unexpected Error.
   Failed to execute (ExtractValuesToTable).
  Failed at Mon Sep 30 19:49:28 2019
  (Elapsed Time: 11 hours 31 minutes
  15 seconds)

This happened numerous times and I don't see what the problem could be, since I performed the same task with much larger files in the past.

I also know that the tool is working at first since I can open the
table and I see entries that look as expected - if I don't open the
table the tool might just crash as well.
I also made sure they are
in the same CRS.
There are only polygons which are completely inside
the raster layer.

So what could be the problem here?


texts:
I'm exporting from PostGis to GeoJson, but appear to have no control over the table it retrieves. Have tried the following methods:

  .\ogr2ogr -f "GeoJSON" "C:\Temp\output.geojson" "PG:host=**** port=**** dbname=**** user=**** password=****" -sql "SELECT * FROM with_id"

SELECT from table with_id failed, no such table/featureclass.

  .\ogr2ogr -f "GeoJSON" "C:\Temp\output.geojson" "PG:host=**** port=**** dbname=**** user=**** password=****" -sql "SELECT * FROM public.with_id"

Unable to open secondary datasource `public' required by JOIN.

  .\ogr2ogr -f "GeoJSON" "C:\Temp\output.geojson" "PG:host=**** port=**** dbname=**** user=**** password=****" with_id

Couldn't fetch requested layer 'with_id'!

  .\ogr2ogr -f "GeoJSON" "C:\Temp\output.geojson" "PG:host=**** port=**** dbname=**** user=**** password=****" public.with_id

Couldn't fetch requested layer 'public.with_id'!


texts: Seeking keyboard shortcut in QGIS for opening attribute table?

texts:
Is there a keyboard shortcut for opening an attribute table for a layer in QGIS ?
I'd like to avoid having to right-click on my layer and then use the drop-down menu to select "Open Attribute Table".


texts: how can openlayers have mobile.css file as style

texts:
I try to make openlayers work in mobile devices, so it can go fullscreen.

I cannot just put width and height to be 100%.

Here is why

 //contains map, forms, images etc

//contains just map's div

//map

//the css
.content{
    padding-left:5%;
    padding-right:5%;
    position: static;
    margin-top:2%;
    margin-bottom:2%;
    width: 50%;
}

#map_element{
width:100%;
height:100%;
z-index:1;
text-align:center;
position:absolute;
}

#mapcont{
width:100%;
height:100%;
z-index:1;
text-align:center;
position:relative;
padding-bottom:60%;
padding-left:0.1%;
padding-right:0.1%;
}

To have a responsive map , its container must have a height. content div cannot have height, because it will not be responsive. This is why I put mapcont.  Map is responsive , but in mobiles is very small (see image) and cannot go fullscreen, because the width and height are already 100%.

So I was wondering how can I use the style.mobile.css as a default css style?

This CSS file is placed in theme/default.

I try to load the CSS like so

<link rel="stylesheet" href="../theme/default/style.mobile.css" type="text/css">

but I see no change on my site.

I found tutorials on how to overwrite the default styles, or how to create new themes. But nothing on how to use the already existing CSS.

EDIT

My apologies if this is still unclear. I was asking for details , so I can implement the rest of the code by my self. Maybe I should show what I am generally trying to do. The whole idea.

This maybe help you understand better. Keep in mind what I said about mobiles, fullscreen and divs and check the following.

User on pc clicks "the map" in the menu and goes here (map, images, forms in the same page)

User on mobile clicks "the map" in the menu and goes here (just the map, fullscreen)

If user wants to see the images and forms clicks "back button" (right corner) and goes here
(just images forms, no map)

If he wants to see the map again, clicks the "back to the map button"

Thanks again for your time, your answers and your patience


texts: Generating Tokens after ArcGIS Server upgrade?

texts:
Our ArcGIS Server was just upgraded from 10.2.x to 10.6 and after this upgrade was done we started to encounter an issue when generating a token to access the ArcGIS services.
The token generation process calls the service from arcgisurl/tokens/generateToken and we pass the parameters as follows:
string url = ConfigurationManager.AppSettings["ArcGIS_Server_URL"] + "/tokens/generateToken/?f=pjson";

            string data = string.Format("username={0}&password={1}&client=requestip&expiration={2}", username, password, ConfigurationManager.AppSettings["TokenExpirationTime"]);

            url = string.Format("{0}&{1}", url, data);

            // Create a new web request
            CustomWebRequest request = new CustomWebRequest(url, "POST");
            //Get Response
            string requestResult = request.GetResponse();

Under arcgis server 10.2.x this was fine, however with 10.6.x  we get the error:
Generate Token request should not contain username and password in the query string

Also I noticed that if I try to access the token generator url from the browser the form comes up, however when I enter the username and password and requestip it returns invalid if done from an external PC. If I try to do the same thing from on the server itself then the token is generated fine, however I don't think that this has any effect as running the application directly on the server still gives the token generation issue.
I also upgraded the proxy to reflect 10.3.x as there was this issue from the query string for that server however this didn't help either.
Has anyone else encountered this issue and if so is there anything that can be done to resolve it?


texts: Is there function to know raster's size (KB) with python?

texts:
I am doing a script with python. There I have many raster in canvas (qgis) which I would like to know size in MB. I am searching function in python that give me this result.

Anyone can help me?


texts: Leaflet custom icons and popups showing GeoJSON properties with external data?

texts:
How can both custom icons and geojson properties be implemented with external geojson in Leaflet?

Using methods described in answers to two other questions, both using inline geoJSON:

Using custom icons for Leaflet Layer groups?

and

Displaying properties of GeoJSON in popup on Leaflet?

I think I have determined that option onEachFeature has to appear before pointToLayer in the L.geoJSON() constructor.

var myceliumTweets = L.layerGroup();
var myLayer = L.geoJSON();
var geojsonFeature1 = 'stream_mycelium.json';

function loadData(data) {
   let xhr = new XMLHttpRequest();
   xhr.open('GET', data);
   xhr.setRequestHeader('Content-Type', 'application/json');
   xhr.responseType = 'json';
   xhr.onload = function() {
       if (xhr.status !== 200) return
           myLayer = L.geoJSON(xhr.response, {
           onEachFeature: function (feature, layer) {
               layer.bindPopup(''+ feature.properties.text +'' + '
' + '
'+
               feature.properties.carmen.latitude + ', ' + feature.properties.carmen.longitude + '
' + feature.properties.profile + '
' + '
'+
               ''+'Carmen' +'' + '
' +
               'Country: ' + feature.properties.carmen.country + '
' +
               'State: ' + feature.properties.carmen.state + '
' +
               'County: ' + feature.properties.carmen.county + '
' +
               'City: ' + feature.properties.carmen.city + '
' +
               'Resolution Method: ' + feature.properties.carmen.resolution_method + '
');
               },
           pointToLayer: function (feature, latlng){
               return L.marker(latlng, {icon:fungi});
               }
           }).addTo(myceliumTweets);
           };
   xhr.send();
};

loadData(geojsonFeature1);

See it working here.


texts: How can you tell if ArcCatalog Geoprocessing is still running?

texts:
Is there some other way other than the logs to check that ongoing geoprocessing is working correctly?

I am running the Buffer tool with dissolve type ALL on ArcCatalog 10.1. I was expecting the the task would take a long time so I left it running and went home for the weekend.  It has now been running for 5 days and I want to find out if it has crashed or not.

Task Manager says the process (RuntimeLocalServer.exe) is still running, and the using between 900 MB and 3.8 GB or memory (it changes over time) and one core of my processor is working flat out, with the others running at about 10%. I have no other applications running.

ArcCatalog displays the hourglass symbol and no error messages in the results pane. I've also checked the log files which just say that the task was started and contain no error messages.

Extra Details

I'm creating a 3m buffer around about 1.5 million polygons in a feature class of a file geodatabase. Many of them are overlapping and so are being dissolved together when they overlap.

I have ArcGIS 10.1 for Desktop on Windows 7 64 Bit.


texts:
I've been trying to figure out how to use ST_Contains to figure out what points are contained within a polygon.
Here is my query.
SELECT "userId", "firstName", "lastName", "phoneNumber", "name", ST_AsGeoJSON("location")::json AS "location", "email" FROM "user_addresses"
INNER JOIN "users" ON "users"."id" = "user_addresses"."userId"
WHERE active AND (
ST_Contains("location",
    ST_GeomFromGeoJSON('{"type":"Polygon","coordinates":[
        [
            [-93.55697, 45.09112],
            [-93.50679, 45.04322],
            [-93.48239, 45.06378],
            [-93.50238, 45.08242],
            [-93.55697, 45.09112],
        ]
            ]}')
    )
);

The query is currently not returning any data,
"location" is stored as SRID=4326;POINT(-93.501611 45.063169)
I'm a bit of a beginner to all of this.


texts: Mosaicing of multitemporal Sentinel-2 images

texts:
My aim is to combine several Sentinel-2 Level2A images of the same location/orbit number but at different dates of a year.

I would like to have different options to create that mosaicing (e.g., if multiple data per pixel is available than create mean/median/std etc., if no data is available than keep Nodata, of only one dataset has data than keep them).

Is there any tool in ENVI which can be used for such mosaicing processes? Alternativly I could also work in ERDAS, SNAP or ESRI ArcMap/ArcPro.


texts: Trying to use the "Iterate Datasets" tool with .dwg files

texts:
I'm trying to expand on a model created by a co-worker that uses the Point and Polyline CAD feature classes within a .dwg to perform a number of geoprocessing tasks.  My goal is to be able to have the location of the .dwg be a parameter for the model and run the model for all of the .dwg files in the folder.  I thought that the "Iterate Datasets" tool would do exactly what I'm trying to do but for some reason the model hits the "Iterate Datasets" tool and then quickly finishes successfully without anything in the rest of the model running.  It acts like it looks in the folder, finds nothing to run, and therefore finishes successfully.  Does this tool not work with .dwg datasets?

The CAD data goes into two separate sub-models after it is iterated.  I rearranged things a little hoping to give a better perspective in the full image.

The only parameter going into the iterator is the folder for the location of the .dwg files.

@Hornbydd, also, this is a bit different now as I'm trying to parse the filename to better name the output folders/gdb/outputs to suit my needs.  Sill having the same problem though.


texts:
I'm trying to have some very basic functionality with vector points with Openlayers 2.12 and I can't figure out how to track a click event on a point.

I create a map using the following code and add one point:

var map;
function init(){
    map = new OpenLayers.Map('map');

    var gsat = new OpenLayers.Layer.Google(
        "Google Satellite",
        {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22}
    );

    map.addLayer(gsat);

    var vectorLayer = new OpenLayers.Layer.Vector("vector");

    vectorLayer.addFeatures(
        new OpenLayers.Feature.Vector(
            new OpenLayers.Geometry.Point(-57.635021, -25.276987).transform(
            new OpenLayers.Projection("EPSG:4326"),
            map.getProjectionObject()
            )
        )
    )

    map.setCenter(new OpenLayers.LonLat(10.2, 48.9).transform(
        new OpenLayers.Projection("EPSG:4326"),
        map.getProjectionObject()
        ), 2
    );

    map.addLayer(vectorLayer);

Then, I've tried to identify when the point is clicked:

    // detect when points are clicked
    pointClick = new OpenLayers.Control.GetFeature(vectorLayer, OpenLayers.Handler.Point,
    {eventListeners:{"click": clickNotice}};
    map.addControl(pointClick);

    // detect point events
    function clickNotice(evt) {
        alert('Point clicked');
        point.deactivate();
    }

}

This is not working, (although I don't get any errors in Firebug), but I'm not sure of the correct syntax for handling points.


texts: Converting SpatialLinesDataFrame to SpatialPointsDataFrame in R?

texts:
In R I have a SpatialLinesDataFrame (for example from a coastline shapefile) and I want to convert it to a SpatialPointsDataFrame. I want to keep all the vertices from the line shapefile to become point shapes in the point shapefile. I also need each point to have the attributes from the original line. For example if the SpatialLinesDataFrame has a line with 5 points and a "name" attribute, then all 5 new points should retain the attribute value.

I figured out some R code to convert the geometry part:

library(rworldmap)
library(sp)
data(coastsCoarse)
lineshape = coastsCoarse

point_coordinates = c()
i = 1
nLines = nrow(lineshape)
for (i in 1: nLines) {
  line1 <- lineshape[i,]@lines[[1]]@Lines[[1]]
  line1coords <- line1@coords
  point_coordinates = rbind(point_coordinates, line1coords)
}

pointshape <- data.frame(x=point_coordinates[,1], y=point_coordinates[,2])
coordinates(pointshape) <- ~x+y

#test result
plot(lineshape)
points(pointshape)

#now how do I transfer the attributes from the lineshape to the pointshape?

I'm stuck on how to transfer the attributes from the SpatialLinesDataFrame to the SpatialPointsDataFrame. I want my script to be reusable for different kinds of SpatialLinesDataFrame variables and so I don't know beforehand the number, names and data types of the attributes.

How can I transfer the attributes?


texts: Filter symbols in QGIS similar to labels for anti-collision?

texts:
I am using the GeoNames data in a PostGIS 2.5 table to populate cities on a map in QGIS, but with millions of entries it is very crowded with overapping making it useless without some sort of filtering. The QGIS layer is points only, and I have various symbols symbols for hamlet, village, etc… classified by population size, and they are all 1mm to 3mm in diameter.

I would like to some how filter the items that get shown in a way analogous to how QGIS draws labels for a layer, with some symbols not getting drawn if they would collide with other symbols. I would also like that to be done with priority for population size rather than the seemingly random way QGIS does it for labels of equal priority in the same layer.
At first I though that scale-based rendering would work, but that means in order to filter out nearby cities like Mississauga or Hoboken so they don't overlap with Toronto and New York at smaller scales, smaller cities like Red Deer or Rapid City would not get shown even though there is plenty of room for their symbols.
I'm imagining the algorithm go something like this:
for each place in list of places in order of most population to least
  get map location where place would be drawn
  have any other symbols already been rendered within a 5mm diameter of that map location?
  if no then
    render this symbol
  else
    skip to next place

In that map portion shown above you can see how there is a chain of six overlapping circles extending out to the left and down from the back dot like a backwards 'J'- large, large, large, small, large, small. Only the large circle between the two small ones would get rendered even though the other large circles indicate cities of larger population since it's the first one to be far enough away from the back dot. If I zoom in to a higher scale, this hypothetical algorithm would let all six circles get rendered once the symbols would no longer be inside each other's 5mm diameter.
I know I can control rendering order in QGIS, but can't seem to figure out the "render or not" part. I found various functions in plugins that do overlaps, but since this is a point layer nothing ever actually overlaps to do it that way. Even if it did, those functions look at all features and I need it too look at only the already rendered stuff to make the decision.
I suspect that this will require Python, but I don't know it that well or what the QGIS internal Python API offers to have an idea of how to go about this. Is this variable rendering even possible in QGIS?


texts:
Attempting to install rgeos package fails.  Detailed error logs as follows:-

> install.packages("rgeos", repos="http://R-Forge.R-project.org")
Installing package(s) into ‘/Users/calvin/R/x86_64-apple-darwin11.2.0-library/2.14’
(as ‘lib’ is unspecified)
Warning: dependency ‘stringr’ is not available
trying URL 'http://R-Forge.R-project.org/src/contrib/rgeos_0.1-15.tar.gz'
Content type 'application/x-gzip' length 166974 bytes (163 Kb)
opened URL
==================================================
downloaded 163 Kb

ERROR: dependency ‘stringr’ is not available for package ‘rgeos’
* removing ‘/Users/calvin/R/x86_64-apple-darwin11.2.0-library/2.14/rgeos’

The downloaded packages are in
    ‘/private/var/folders/4l/rd4wl50x52n504t2vzqskhhc0000gn/T/RtmppBNwSt/downloaded_packages’
Warning message:
In install.packages("rgeos", repos = "http://R-Forge.R-project.org") :
  installation of package ‘rgeos’ had non-zero exit status
> install.packages("stringr", repos="http://R-Forge.R-project.org")
Installing package(s) into ‘/Users/calvin/R/x86_64-apple-darwin11.2.0-library/2.14’
(as ‘lib’ is unspecified)
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
  package ‘stringr’ is not available (for R version 2.14.0)

Suggestions?


texts: Delete content from vector files stored in Geopackage

texts:
I have packaged a group of vector files (shapefiles, 30-ish) into a Geopackage database including styles, with a view to being able to use that GPKG as a standardised template with containers for future projects - to standardise layer names and default styles etc.

The individual layers still contain the entities associated with the last project from which they are archived.  The layers are a mix of point, line and polygon formats.  Is there a simple way to delete all entities within each layer within the GPKG, or do I need to open the layers within QGIS, edit (delete) and save them back to the GPKG individually?


texts: How to use IDatabaseConnectionInfo3 Interface ArcObjects!

texts:
I try to retrive the connected users from an entreprise geodatabase using ArcObjects .

first I tried to cast My workspace to IDatabaseConnectionInfo3 but it throws

  ClassCastException

IWorkspace workspace = sdeWorkspaceInstance();
IDatabaseConnectionInfo3 databaeseConnectinfo = (IDatabaseConnectionInfo3) workspace;
IEnumUserInfo listUserInfo = databaeseConnectinfo.getConnectedUsers();

Any idea how to use the IDatabaseConnectionInfo3 interface ?


texts: How do I get the "clip with another polygon layer" tool/plugin from QGIS version 2.14.1

texts: Hello I'm looking for help finding and downloading a plugin/tool from a previous version of QGIS (2.14.1). It's called the "clip with polygon from another layer" and I would like to download it for the latest version of QGIS (2.18.14). Whenever I've searched plugins on the QGIS desktop it doesn't show up and I've also looked on the QGIS plugin website online and it isn't in their either to my knowledge. If anyone has any idea of how i can get this plugin/tool on the latest version it would be very helpful.


texts: suitable column Data-type for geojson data storage and processign in PostGis


texts: Get Radius in epsg 4326 with Openlayers 3

texts:
What is the best method for getting the radius of a circle in Openlayers 3 in epsg 4326 (geodesic)?  I'm using this example for a base and trying to draw a circle while getting the radius, but using proj 4326 is proving difficult.
Here's my code:

var pointerMoveHandler = function(evt) {
   if (evt.dragging) {
     return;
  }
var tooltipCoord = evt.coordinate;

  if (sketch) {
    var output;
    var geom = (sketch.getGeometry());
    if (geom instanceof ol.geom.Cirlce) {
      output = formatCircle(/** @type {ol.geom.Circle} */ (geom));
      tooltipCoord = geom.getCenter().getRadius();

    measureTooltipElement.innerHTML = output;
    measureTooltip.setPosition(tooltipCoord);
  }
};

var sourceProj = map.getView().getProjection();
    var geom = /** @type {ol.geom.Circle} */(Circle.clone().transform(
        sourceProj, 'EPSG:4326'));
    var center= geom.getRadius;
    radius = Math.abs(wgs84Sphere.geodesicArea(center));
var output;
  if (area > 10000) {
    output = (Math.round(area / 1000000 * 100) / 100) +
        ' ' + 'km2';
  } else {
    output = (Math.round(area * 100) / 100) +
        ' ' + 'm2';
  }
  return output;
};

radius = Math.abs(wgs84Sphere.geodesicArea(center)); isn't right, but I'm not sure what to put there....


texts: Changing the symbology of a .geojson file in QGIS, through the Python console

texts:
I have tried to look through the PyQGIS developer cookbook and change the symbology of a point file to a black circle but unfortunately every time I add it, the symbology does not change. It should just be a simple process but it does not seem to work. The first line of code I have written so far adds the .geojson file to the project and the second one aims to change the symbology:
vlayer = iface.addVectorLayer("path_to_.geojsonfile", "Earthquake Event: ", "ogr")
if not vlayer:
    print("Layer failed to load!")

symbol = QgsMarkerSymbol.createSimple({'name': 'circle', 'color': 'black'})
layer.renderer().setSymbol(symbol)
# show the change
layer.triggerRepaint()


texts: How to pass a variable (FilesToAdd) from VBA script (interface) to a Mapbasic program?

texts:
Option Explicit

Public FilesToAdd As Integer

Private Sub TEST_Click()

FilesToAdd = 1

Dim b As Variant
Dim Prog As String
Dim File As String

Prog = "C:\Program Files\Mapinfo\Professional\Mapinfow.exe"
File = "C:\Mapinfo\Combine\SelectCarrier.MBX"

b = Shell(Path + " " + File, vbNormalFocus)

End Sub


texts: Changing metadata using ArcPy of ArcGIS Pro?

texts: Is there any ArcPy function or classes which allow us to change map and layout metadata within ArcGIS Pro?


texts: I'm trying to import an SHP file containing the official map of rural bus routes in Ireland (https://www.transportforireland.ie/transitData/Rural_Transport_Services.zip). The file loads without errors but instead of expected lines qGIS plots strange polygons (see image) which by no means represent public transport routes. What am I missing?


texts: Create shortest route for a set number of units

texts:
I am running a field survey and would like to create routes for the surveyors. I have the building lots and road network for a particular city (road in red, lots in blue, with a Google Maps basemap):

Each of the lots may have more than one apartment unit, and the number of units is tracked as an additional column of the attribute table.

For routing purposes, I would like to generate routes of 50 units each for the surveyors. Is there a way to create a script (preferably in ArcGIS) that can iterate through the entire lot shapefile to create the most compact routes of 50 units at a time. Ideally, the output would be a new column in the lot shapefile that designates for each lot the route that it is assigned to.


texts:
I am trying to create a public WFS service with GeoServer from a PostGIS SQL view, but this does not load in QGIS. If I publish from a PostGIS table the WFS service is loaded normally.

Is it possible to publish a WFS service from an SQL view?


texts:
There are quite some nice possibilities to analyse spatial data in R and in the context of my current project I would like to use R more frequently to do that.

Until now I plot my maps with ggplot2 package that brings a lot of practical tools to plot and explore data. Still my maps are not nearly as good looking as the ones I used to plot in ArcGIS.

So for publishing purposes I wonder if there are any good tutorials, books, practical tips, packages etc. that I could use to experiment a little bit and make my maps more good looking?


texts: AssetFinda Plugin projection issues in QGIS

texts:
I am currently trying to overlay raster and vector data in a common projection as the Asset data which is using NZGD:2000 EPSG:2193. I have an ecw aerial image that initially loads as a generated CRS. There are also parcels and road vector data. I reprojected the vector layers using the 'save as' method and I reporjected the ecw image using the wrap tool. They I then load them into a new map project. I set the CRS for the map project to be epsg:2193 I also enable on the fly project. After that I load up the AssetFinda layers(which are sourced from an Sqllite.db generated from view tables from a SQL database) and unfortunately they still don't overlay with the other layers.

I had a thought that perhaps they are in the same projection but have different extents. However, they should have the same extent. Basically I have tried all the standard reprojection techniques and I was hoping someone here can think of something I may have overlooked.


texts: Split list by chunks Earth Engine with mean

texts:
I am calculating mean temperature data based on historical data. I need to calculate the mean through the years (jan-jan-jan-jan...) and the means within a year (jan-feb-mar-may-jun...).

This is what I have but it is not working at all. I am new coding in jave and Google Earth Engine.

How can I approach to solve this issue?

// Load TerraClimate data
var dataset = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE').select('tmmx');

var maximumTemperature = dataset
    // Define reference conditions from the first 25 years of data.
                  .filter(ee.Filter.date('1990-01-01', '2015-12-31'));
   // Sort chronologically in descending order.
                  //.sort('system:time_start', false);

print('maximumTemperature', maximumTemperature);

// list for filter iteration
var months = ee.List.sequence(1,12);


texts:
I am exporting a single image and a collection of Landsat raw images from Google Earth Engine using Python, but all my images are in black and white. Should they be like that? Isn't there a way of getting coloured Landsat images. I am very new.
Below is my code:

# Import libraries.
import ee
import webbrowser

# Trigger the authentication flow.
# ee.Authenticate()

# Initialize the library.
ee.Initialize()

# Access a collection
collection = 'LANDSAT/LC08/C01/T1' #Landsat 8 Collection

#Landsat_composite in Nithi area
area_nithi = list([(37.288, 0.172), (38.333, 0.172), (38.333, -0.568), (37.288, -0.568), (37.288, 0.172)])
area_nithi = ee.Geometry.Polygon(area_nithi)
time_range_nithi = ['2017-01-01', '2017-12-31']

collection_nithi = ('LANDSAT/LC08/C01/T1')
print(type(area_nithi))

#Methods from climada.util.earth_engine module
def obtain_image_landsat_composite(collection, time_range, area):
    """ Selection of Landsat cloud-free composites in the Earth Engine library
    See also: https://developers.google.com/earth-engine/landsat

    Parameters:
        collection (): name of the collection
        time_range (['YYYY-MT-DY','YYYY-MT-DY']): must be inside the available data
        area (ee.geometry.Geometry): area of interest

    Returns:
        image_composite (ee.image.Image)
     """
    collection = ee.ImageCollection(collection)

    ## Filter by time range and location
    collection_time = collection.filterDate(time_range[0], time_range[1])
    image_area = collection_time.filterBounds(area)
    image_composite = ee.Algorithms.Landsat.simpleComposite(image_area, 75, 3)
    return image_composite

#Application to examples
composite_nithi = obtain_image_landsat_composite(collection_nithi, time_range_nithi, area_nithi)

#Selection of specific bands from an image
nithi_band = composite_nithi.select(['B4', 'B3', 'B2'])

nithi_bandVis = {
  min: 0.0,
  max: 30000.0,
}
print(composite_nithi.getInfo())

print(type(nithi_band))

print(type(nithi_bandVis))

def get_region(geom):
    """Get the region of a given geometry, needed for exporting tasks.

    Parameters:
        geom (ee.Geometry, ee.Feature, ee.Image): region of interest

    Returns:
        region (list)
    """
    if isinstance(geom, ee.Geometry):
        region = geom.getInfo()["coordinates"]
    elif isinstance(geom, ee.Feature, ee.Image):
        region = geom.geometry().getInfo()["coordinates"]
    elif isinstance(geom, list):
        condition = all([isinstance(item) == list for item in geom])
        if condition:
            region = geom
    return region

region_nithi = get_region(area_nithi)

def get_url(name, image, imageVis, scale, region):
    """It will open and download automatically a zip folder containing Geotiff data of 'image'.
    If additional parameters are needed, see also:
    https://github.com/google/earthengine-api/blob/master/python/ee/image.py

    Parameters:
        name (str): name of the created folder
        image (ee.image.Image): image to export
        scale (int): resolution of export in meters (e.g: 30 for Landsat)
        region (list): region of interest

    Returns:
        path (str)
     """
    path = image.getDownloadURL({
        'name':(name),
        'scale': scale,
        'region':(region)
        })

    webbrowser.open_new_tab(path)
    return path

#url_nithi = get_url('nithi', median_swiss, 900, region_swiss)
url_nithi = get_url('nithi0.1', composite_nithi, 900, region_nithi)
#url_landcover = get_url('landcover_swiss', landCover_layer, 100, region_swiss)

#For the example of nithi, due to size, it doesn't work on Jupyter Notebook but it works on Python
url_nithi2 = get_url('nithi0.2', nithi_band, nithi_bandVis, 1000, region_nithi)

print(url_nithi)
print(url_nithi2)


texts: Using the expression builder in the Atlas function in QGIS to group features?

texts:
I am trying to show only record in the atlas for the following attribute table...

ogc_fid | id |     phasename       | exampleinfo
----------------------------------------------------------
   2    | 0  | 1- 3 Lower Hillgate | example information 1
   2    | 0  | 1- 3 Lower Hillgate | example information 2
   2    | 0  | 1- 3 Lower Hillgate | example information 3
   2    | 0  | 1- 3 Lower Hillgate | example information 4
   2    | 0  | 5 Lower Hillgate    | example information 1
   2    | 0  | 5 Lower Hillgate    | example information 2
   2    | 0  | 5 Lower Hillgate    | example information 3
   2    | 0  | 3 Lower Hillgate    | example information 1

So in my atlas I only want their to be 3 records.

Is their a simpler way to do this as opposed to making a copy of the data and having to then use the combined data to make an atlas?

The reason why I am doing this is because the example information is different and i want to show the attribute table data associated to the phasenamein an attribute table shown below the atlas map window I am trying to create.


texts: Drawing tools in ArcGis JS not always registering clicks in IE 11

texts:
When constructing a toolbar such as in the example given at http://developers.arcgis.com/javascript/samples/toolbar_draw/ there is a slight issue when you're trying to place to first point of any drawing tool. If you are too fast (it seems to be about 2 seconds) the click won't be registered at all. Has anyone noticed this, and found a work-around or something? There is another example here http://developers.arcgis.com/javascript/samples/graphics_add/ where this problem doesn't exist, but I haven't been able to figure out what's causing this bug.

This also goes away if you set the document mode to IE 9 and also it doesn't exist at all in Chrome.


texts: Cache generation failed. Progress gets stuck at 0%

texts:
I'm trying to publish a map service containing a large .tif raster (48GB).

I managed to publish a small cached raster successfully (for testing), but when I tried to publish the large one, the cache generation got stuck at 0% even though it kept saying it's in progress, and again at 24%. Now I can't even publish the small one because the application is not responding right after I select "Publish".

Details:

It's an ArcGIS Server Enterprise 10.1 on a Windows Server 2008 OS machine.

the machine is on a Microsoft Azzure cloud.
the machine has 4 CPU cores (so I set the number of instances to 5).
it's a raster TIF file with LZW compression.
I chose 5 scales, and tried different formats for the tiles.

Thanks for your help


texts:
I am not very computer literate; I get by, but barely. I have been using google earth for several years to plot ancient river channels and old mining claims from MRDS on Google Earth. I have (used to have) pretty extensive maps drawn from all sorts of sources- until google earth crashed, and I lost several hundreds of hours of work- in several instances (I am a slow learner about GE). So I "get" GE (except for the crashing part) and I had multiple layers on it until it ALL went away.

I looked around and it sounds like QGIS may work better, but I have spent a lot of hours trying to figure out how to get QGIS started- despite spending hours looking at "QGIS 101" videos on Youtube. I am very frustrated as this is a GIS terminology is a foreign language to me and apparently the conventions and terminology are foreign to me.

What I am looking to do is to simply add several layers onto either a topo and a satellite image. That's it. It would be nice if the program was set up simply like GE, but would have more defined layers (which I am sure EVERY program has).

I have spent hours trying to get appropriate maps loaded and I am being confounded with too many options and terminology in the windows that I don't understand. I have tried to follow several videos online, and the screens shown on the videos are different from the ones I am viewing. I feel I am spending way too much time, just figuring out how to set up a simple program.

Rather than waste hours of my time trying unsuccessfully to set this up, is there any "service", where a person can be "coached" on how to set up a program like QGIS? My needs are simple, and I assume as soon as I get the program set up, I can start using it, but I am having issues setting it up, despite reading and viewing tons of stuff.

I am stuck on how to download topo maps and how to know about and to get the plugins that I need to get started. It would be nice to find someone to just walk me through the set up with what I need to get started, if QGIS is my answer. Is there anywhere I can go to get this help?


texts: Using Copy / Paste when editing attribute tables in QGIS

texts:
I have a number of polygon attribute tables with columns that I want to update with new information.

Is there a way of (relatively) quickly editing info in attribute tables in QGIS?

For example, such as using Control C and V such as in MapInfo. Currently seem to need to double click/right click/copy, move to cell, double click/right click/paste each individual cell I want to change - would like something simple rather than going down the field calculator, filter, etc route. With points no problem; I usually throw them out to Excel, edit and then bring back in again.


texts: OL3: Is there a getUrl function for ol.source.XYZ?

texts:
Problem: I am designing an ol3 plugin that needs to copy the URL parameter from ol.source.XYZ/ol.source.TileWMS to ol.layer.Tile on map initialization.

TileWMS: As per the docs, this is easily done with the following.
Note: not actual code, simplified for example. Assumes all layers are TileWMS with no layer groups.

var layers = this.getLayers().getArray();
for (var i = 0; i < layers.length; i++) {
  layer = layers[i];
  layer.set('url', layer.getSource().getUrls()[0]);
}

XYZ: In the docs, I see a setter but no getter. Currently I'm using layer.getSource().setTileLoadFunction() to gain access to the imgSrc parameter, which is the actual URL orf a loaded tile, and genericize it back into a template URL. This works with a single layer paired with one source, but breaks when multiple layers share the same source.
Note: Same assumptions as above.

var layers = this.getLayers().getArray();
for (var i = 0; i < layers.length; i++) {
  layer = layers[i];
  var oldGetTileLoadFunction = layer.getSource().getTileLoadFunction();
  layer.getSource().setTileLoadFunction(function (imgTile, imgSrc) {
    imgTile.getImage().src = imgSrc;

    layer.set('url', MyApp.utils.generateXYZUrl(imgSrc));
    layer.getSource().setTileLoadFunction(oldGetTileLoadFunction);
  });
}

Is there any way to have something like layer.getSource().getUrl() for XYZ?
Note: must use native javascript for solution, no external libraries allowed


texts:
I've been adding easements, which typically have a fixed width. Are there any existing tools or scripts that would allow me to define one edge and automatically create the other? The buffer tool is poor, the line offset tool is for lines.

I want to be able to use the same tools (sketch, intersect, endpoint arc) and have actual curves,not point approximations (what buffer does).


texts: Getting clipped Landsat RGB image to render correctly in ArcMap

texts:
I downloaded Landsat 8 tiff bands 1-7 into Arcmap 10.7 and then used the Composite Band tool to create a RGB raster composite image. In order to have the Landsat image render correctly I had to reassign the bands so that Red was band 4, Green was 3, and Blue was 2.

I then used the Clip tool to clip the raster using a polygon layer, but the resulting RGB raster is rendering incorrectly even though the bands are assigned the same way. The clipped raster appears off-color and tinted in regions, and is not the same image quality as the first raster so I think it has to do with the RBG band assignment.

But I've tried different band assignment combos thinking that the bands were shuffled during the extraction and band 3 isn't really 3, etc., but no luck. Is there a solution to render the clipped raster like the original, or would it work if I re-assign bands AFTER I clip to raster for it to render properly?


texts:
I wanted to plot some raster-data with mapview. I always get this greyish bounding box around my data though. Anybody an idea of what this might be?
I just call mapview(<raster>)


texts: accent-insenstitive text in CartoDB - PostGreSQL unaccent function missing?

texts:
I am trying to do accent-insensitive text matching, but PostGreSQL's unaccent function doesn’t appear to be available

SELECT unaccent('Hôtel');

function unaccent(unknown) does not exist

Is this function supported somehow in CartoDB? If not, is there an alternative?


texts: Python fastest way to load overlapping part of 2 rasters with different origin to np.array

texts:
I got a problem with gdal.wrap - it works REALLY slow. On 10m raster of Poland, matching the origin of 2 rasters takes like 30 minutes.
Is there an easy and fast way to take 2 rasters (of different origins) and to load their overlapping part into the numpy arrays?
Rasters are having the same resolution and crs.
    def run_mosaic_script(self, resample_alg: any = None, x_res=13.9, y_res=13.9, compress_algorithm="LZW",
                          dstSRS="EPSG:32634") -> bool:
        """
        Method for run mosaic script

        :param compress_algorithm: compression algorithm example LZW
        :type compress_algorithm: str or None
        :param resample_alg: resample algorithm
        :type resample_alg: any
        :param x_res: x resolution of pixel
        :type x_res: numeric (int or float) or None
        :param y_res: resolution of pixel
        :type y_res:numeric (int or float) or None
        :return: True if correct, False if incorrect
        :rtype: bool
        """
        creation_option: List[str] = ["BIGTIFF=YES"]
        if compress_algorithm is not None:
            creation_option.append("COMPRESS=" + compress_algorithm)
        gdal_warp_options = gdal.WarpOptions(srcNodata=0, dstNodata=0, multithread=True,
                                             warpOptions="NUM_THREADS=ALL_CPUS", xRes=x_res, yRes=y_res,
                                             dstSRS=dstSRS,
                                             resampleAlg=resample_alg, creationOptions=creation_option,
                                             cutlineDSName=self.shp_file.__str__(), cropToCutline=True,
                                             overviewLevel='AUTO')
        image_to_mosaic: list = list()
        for file_name in self.input_image_list:
            image_to_mosaic.append(file_name.__str__())
        gdal.Warp(destNameOrDestDS=self.output_image.__str__(),
                  srcDSOrSrcDSTab=image_to_mosaic,
                  options=gdal_warp_options)
        return True

I am just cropping it to some simple extent. Both of the rasters.


texts: I am trying to run the Grouping Analysis tool from Esri's ArcGIS Spatial Statistics/Mapping Clusters toolbox.  I have run the tool successfully on the same input layer using a different spatial weights matrix, but when I try to run the tool with an IDW spatial weights matrix the tool fails.  I have not seen this type of error before.  See image...


texts: Writing shapefile of multiple geometry using Python

texts: How to make a shapefile containing points, polygon, and line using python shapely or pyshp?


texts:
I would like to add dimensions to my QGIS map. It can be in the editing/drawing space or the print composer. I would like to see how long a line is, for example 3,5 m.

I tried a lot of different things, but I can't figure out how to add dimensions to my map.

Is there a way to add dimensions to a QGIS map?


texts: Calculate from tile URL path a TileCache URL path

texts: I am searching for a way how to calculate in PHP (there for OpenLayers are not suitable for me) from tile path /zoom/x/y.png a Tilecache URL like 14/000/009/107/000/010/749.png


texts:
I am trying to create a simple webmap and I want to add the layers of a map that I published to the GIS server. But the rest URL from the published map do not display these layers. Anyone know what the problem is?

I noticed that I was missing the section "Supported Interfaces:   REST   SOAP   WMS" as you see from the picture, could that be the reason?

I am very new to this.

Below is a snipped of code that I used to test the rest url.

  require([
    "esri/map",
    "esri/layers/ArcGISDynamicMapServiceLayer",

  ], function (
    Map, ArcGISDynamicMapServiceLayer) {

    map = new Map("mapDiv", {
    basemap: "topo",

      sliderOrientation : "horizontal"
    });

    var dynamicMapServiceLayer = new ArcGISDynamicMapServiceLayer("rest url link");

    map.addLayer(dynamicMapServiceLayer);


texts: What files do you need to begin the process of creating a public transportation multimodal network in QGIS?

texts:
I have never created a multimodal network from scratch before. I want to use General Transit Feed Specification (GTFS).txt files for a New York MTA public transportation network.

What files do I need to begin the process before I even ask about the steps to be performed in QGIS? I already have a lions street file.


texts: Batch export of image collection giving consecutive numbers in filename

texts:
I want to use consecutive numbers in the name of the exported images from a image collection. I want to filter the MODIS image collection within a specific date range and region for the worst EVI images, rank them accordingly to their statistics and export the worst 5 images with their position in their filename. So the image with the worst EVI-mean value gets the filename Top 1 worst EVI, the second worst Top 2 worst EVI and so on. So far, I am just able to give all the 5 worst images the name Top 5 worst EVI, but I would like to make that more specific.

var ModisFiltered = MODIS.filterDate('2019-01-01', '2020-05-22').filterBounds(ROI)
  .map(function(image){return image.clip(ROI)});

var maxPixels = 1e9; //Maximum amount of pixels for ee.Reducer calculations
var amountWorstYears = 5;

//Calculate Mean of EVI
var FinalCollection = ModisFiltered.map(function(toSVImean){
  var selectSVIband = toSVImean.select(['EVI']);
  var calculateEVImean = ee.Number(selectSVIband.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: ROI,
    scale: 1000,
    bestEffort: true,
    maxPixels: maxPixels
  }).values().get(0));

  var addSVI_mean = selectSVIband.setMulti({
    EVI_Mean: calculateEVImean
  })

  addSVI_mean=ee.Image(addSVI_mean)
return addSVI_mean.addBands(toSVImean);
})
print ('Image Collection with EVI Mean:',FinalCollection)

//Sort and Filter for lowest SVI
var select = FinalCollection.select(['EVI']).limit(amountWorstYears, 'EVI_Mean');

//Add amount of requested "worst SVI products" to the properties if each image
var selectionWithNumbers = select.map(function(addAmount){
  return addAmount.set({Position_worst_EVI: amountWorstYears})
});
print('Selection with worst EVIs',selectionWithNumbers)
Map.addLayer(selectionWithNumbers.first());

//Calculate consecutive position numbers (suggestion to provide consecutive numbers)
function CalcRank (current, previous){
  var previousNumber = ee.Number(previous//.getNumber('Position_worst_EVI')
  )
  return current.set({Position_worst_EVI: previousNumber})
  return previousNumber.add(1)
}
var rankedSviImages = ee.ImageCollection(selectionWithNumbers.iterate(CalcRank, 0))
print(rankedSviImages)

//Export Image Collection
var batch = require('users/fitoprincipe/geetools:batch');

batch.Download.ImageCollection.toDrive(selectionWithNumbers, "SVI_ImageCollection",{
  name: 'Top{Position_worst_EVI}_worst_SVI_{system_date}',
  scale: 1000,
  region: ROI
});

I am trying to solve my problem with an iteration, but as this is something new for me, I do not know how to apply that iteration correctly.

//Calculate consecutive position numbers (suggestion to provide consecutive numbers)
function CalcRank (current, previous){
  var previousNumber = ee.Number(previous//.getNumber('Position_worst_EVI')
  )
  return current.set({Position_worst_EVI: previousNumber})
  return previousNumber.add(1)
}
var rankedSviImages = ee.ImageCollection(selectionWithNumbers.iterate(CalcRank, 0))

Here is the GEE Code


texts: Automating Print layout with export map canvas as image using PyQGIS

texts:
I would like to automate two steps in my QGIS Print Layout by using PyQGIS
I have my print layout already created, and I was trying to open it by this code:
composerTitle = 'Strand3'

project = QgsProject.instance()
projectLayoutManager = project.layoutManager()
layout = projectLayoutManager.layoutByName(composerTitle)

which comes from this query:
How to select a specific print layout (by name) in QGIS 3?
The code didn't work as I wanted. I just wanted to load my print layout called "Strand3" from the layout manager. However is not as important as the stuff below, because I can load it once manually and keep it open.
Next, I would like to do an automatic Set Map Extent to Match Main Canvas Extent.
I found some solution here:

Setting the extent for a QGIS map using python
Set print composer to map canvas extent using python
Map extent matching permanently main canvas extent

but it doesn't work, as they propose the definition of map rectangle. I simply need my map refreshed, what was considered here:

QGIS print layout - map refresh not working

I would like to have this step done automatically as well as the last one - exporting map to PNG.
Since I've managed with setting the fixed DPI resolution for the map, I want to skip the window below:

when possible.
I found the code like this:
 .exportToImage(fn, QgsLayoutExporter.ImageExportSettings())
 exporter.exportToPdf(fn, QgsLayoutExporter.PdfExportSettings())

and solution like this:

https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/composer.html#output-using-print-layout

but I don't know how to implement Qgis Layout Exporter.
Anyway my final code looks like this:
 comp = iface.createNewComposer()
 comp.composition().loadFromTemplate(myDocument)

 map_item = comp.composition().getComposerItemById('Strand3')
 map_item.setMapCanvas(canvas)
 map_item.zoomToExtent(canvas.extent())

 base_path = os.path.join(QgsProject.instance().homePath())
 png_path = os.path.join(base_path, ".png")

 exporter = QgsLayoutExporter(layout)
 exporter.exportToPdf(png_path, QgsLayoutExporter.ImageExportSettings())

and I am getting an error: 'QgisInterface' object has no attribute 'createNewComposer'.
How can I automate these 2 things in my QGIS print layout?


texts: Identify polygons intersected by lines in a different shapefile

texts:
I have a shapefile with lines and a different shapefile with polygons. I want to identify the id (attribute) of the polygons intersected or overlapped by each line, to then make some additional calculations.

I was thinking about finding the number of intersections but that wouldn't give me the information regarding the polygon id.

I think using the intersect tool from QGIS with the id as selected feature could work, but I can't get it working.

I'm using QGIS 2.18.12


texts: Cross sections of two surfaces in Civil 3D

texts:
In Civil 3D I have calculated the volume between two surfaces. I would like to have corresponding cross sections of the first and second surface for all the study area.
"Create superimposed profile" cannot solve the problem completely since it creates only profiles along one alignment and doesn't cover entire study area.
I have did it in Land but I could not find the tools in Civil 3D.

Any suggestion.


texts: Wrong $area from field calculator (QGIS)

texts:
The area of the polygons in a layer, calculated by the field calculator with "$area" is wrong (much too big) for some polygons. If I draw the concerned polygon again, the area is calculated correctly. It looks like the polygons really have a wrong size, but the centroids are calculated correctly.

edit: Thanks for your replies so far. Here is an example with some polygons. For some, the area is calculated right, for some not (as I said much too big). For example the polygon with the ID 669. If I redraw it, the area is calculated right. The problem is, I have >1400 polygons, and I don't have the time to check and redraw every single one.

Maybe someone got an explanation/solution?


texts:
I am trying to add a button to an ArcPad form that sets the values to Cereal and Corn. This is as far as I have gotten and I am stuck.

I Made apl that has a form. I get the values into the form using the arcpad data manager from ArcMap. it populates the form giving me combo boxes. What I need is a button that will make it easier to choose specific values from the combo boxes. So instead of scrolling down the combo boxes to make my choice and hitting OK. I want to just hit the corn button and then that updates the combo boxes and adds the form and closes it.
In the image you can understand what it is I need to do. This is for field work and the easier it is to insert some crop types the faster it gets done.

'once called updates the fields to add in corn to the form then closes the form
Sub InsertCorn(ByVal layerName, ByVal pageName)

   Dim objEditForm
   Dim objTheLayer
   Dim objEFPageControls
   Set objTheLayer = Application.Map.Layers.item(layerName)
   Set objEditForm = objTheLayer.Forms(1)
   Set objEFPageControls= objEditForm.Pages(pageName).Controls

   objEFPageControls("cboAGRLU_ClassGRP").value = "Cereal"
   objEFPageControls("cboCropType_Varietal").value = "Corn"

   Set objEFPageControls= Nothing
   Set objEditForm = Nothing
   Set objTheLayer = Nothing
End sub

then in the onclick event on the button:

  Call InsertCorn("Corn", "PAGE1")


texts: PyQGIS XYZ Tiles not rendered if jpeg tiles are returned

texts:
In a stand-alone/external PyQGIS app that I'm developing I've noticed that wms/xyz tiles that are returned as jpegs don't render on the map canvas but those that return pngs do.

So for example the Google Maps, layer=m which returns pngs does render fine.

provider='type=xyz&url=http://mt1.google.com/vt/lyrs%3Dm%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=19&zmin=0'

But the Google Satellite, layer=s which returns jpegs doesn't render, nor does the Bing layer which also returns jpegs render.

provider='type=xyz&url=http://mt1.google.com/vt/lyrs%3Ds%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=19&zmin=0'

I've tried adding an image format parameter, format=image/jpeg but it doesn't make any difference.

provider='format=image/jpeg&type=xyz&url=http://mt1.google.com/vt/lyrs%3Ds%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=19&zmin=0'

I used Fiddler on Windows to confirm that there were no network issues with the Google Satellite requests and I can see the jpeg tiles being returned successfully. It's actually how I first noticed that the working layers were returning pngs and the non-working layers were all returning jpegs.

Is there some parameter I haven't set for the raster layer to ensure that jpeg tiles are rendered?

The same wms/xyz tile URI, and without the format parameter works fine in QGIS itself, just not in my stand-alone PyQGIS app.

I'm using QGIS 3.4.

Here is a minimal PyQGIS app demonstrating the issue:

from qgis.core import QgsApplication, QgsProject, QgsRasterLayer
from qgis.gui import QgsMapCanvas, QgsLayerTreeMapCanvasBridge

from PyQt5.QtWidgets import QMainWindow
import os

GOOGLE_MAPS_SOURCE = "type=xyz&url=https://mt1.google.com/vt?lyrs%3Dm%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=19&zmin=0"
GOOGLE_SATELLITE_SOURCE = "type=xyz&url=https://mt1.google.com/vt?lyrs%3Ds%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=19&zmin=0"

class WMSTestApp(QMainWindow):

    def __init__(self):
        super().__init__()

        self.canvas = QgsMapCanvas(self)
        self.bridge = QgsLayerTreeMapCanvasBridge(QgsProject.instance().layerTreeRoot(), self.canvas)

        layer = QgsRasterLayer(GOOGLE_SATELLITE_SOURCE, "Test WMS Layer", "wms")
        print("Raster layer is valid?:", layer.isValid())

        QgsProject.instance().addMapLayer(layer)
        self.canvas.setLayers([layer])

        self.setCentralWidget(self.canvas)

QgsApplication.setPrefixPath(os.environ.get("QGIS_PREFIX_PATH", "C:\Program Files\QGIS 3.4\apps\qgis"), True)
qgs = QgsApplication([], True)
qgs.initQgis()

app = WMSTestApp()
app.show()
qgs.exec_()
qgs.exitQgis()

Fiddler showing the successful transmission of the jpeg tile.


texts: Differences between Geoserver 2.15.0 and 2.11.0

texts:
I was doing a GetFeature request call to a WFS service in Geoserver in json format with my application. I had initially installed 2.11.0, then when I moved to 2.15.0 my app stopped running,

I have discovered that the header response in the json is different depending on the Geoserver version I'm using,

In 2.11.0:
{"type":"FeatureCollection","totalFeatures":56,"features":[

In 2.15.0:
{"type":"FeatureCollection","features":[

So the number of total features is missing,

What Geoserver configuration step am I missing in order to get the same response?


texts: How do I select all overlapping features (in the same layer)?

texts:
I recently moved to QGIS 2.18.3

On my old QGIS when I clicked on a feature with another feature(s) beneath it, in the same layer, it would highlight them all and all beneath my mouse click or box would be selected in the attribute table.

Now on QGIS it only selects the top feature, which restricts how I use QGIS.
Is there a way to select all under a point?
Or at the very least, is the a keyboard button I can press to cycle through the features when I click, like I used to be able to do in mapinfo. - I've tried a few to no avail.


texts: Adding rows of data to a non-spatial table in ArcMap with a SDE SQL Server for Metadata purposes

texts:
I'm attempting to use several non-spatial tables to better organize the large set of feature classes I continually add to a SDE Geodatabase (SQL Server Express).  So I have several tables setup and the relationships between them setup.  Basically the plan is to use these tables to simply add metadata so to speak so the feature classes can be easily reported on.  The tables will show who the owner of the geographic area the data comes from and what type of data it is.  None of this needs to spatial.  Basically I want organization outside of needing to add this metadata type data in every record of the feature class.

However, after creating these non-spatial tables there doesn't seem to be way to add any data to these inside of ArcMap.  The Editor doesn't work since it since there are no spatial layers present.  I've tried arcpy and I only get errors.  Is the only way to use these non-spatial tables is to create some interface outside of ArcMap?  I was think I could use some 3rd Python module to add data directly to SQL Server but have no idea where to find that. Also, how would I create a UI for the Python Scripts?  Any ideas or direction would be appreciated.


texts:
Using esri javascript api version 3.8, I am trying to use some cached tiles from our server at http://hexe.er.usgs.gov/ifhp/will/tiles/ using WebTiledLayer. I believe that our tiles are set up as ${level}/${row}/${col}.png , however my tiles still aren't showing up. Does this have something to do with setting up an extent parameter? I'd like the tiles to show up with the extent of 'xmin': -88.243751, 'ymin': 41.49263, 'xmax': -88.131789, 'ymax': 41.727235

This sample is the only example I've found that uses a web tile cache and it doesn't use the extent parameter but its a world map whereas mine is a small county in Illinois. I haven't found any examples of using the extent option so maybe I just don't have it set up correctly.

Here's my code snippet with the attempt at the extent parameter.

var tiles = new WebTiledLayer("http://hexe.er.usgs.gov/ifhp/will/tiles/overlay10/${level}/${row}/${col}.png", {
    "Extent": "'xmin': -88.243751, 'ymin': 41.49263, 'xmax': -88.131789, 'ymax': 41.727235"
});
map.addLayer(tiles);


texts: From osm.pbf data to mbtiles (pbf vector tiles inside a sqlite3 container)?

texts:
I need to experiment with vector tiles based on OSM data.
To do that, I've downloaded an OSM *.osm.pbf data file from geofabrik.de.
I want now to store these data either as the standard /tileset/{Z}/{X}/{Y}.vector.pbf files, or as a single sqlite3-flavoured *.mbtiles file (better).
Do you know a straightforward way to proceed using command line tool(s) on Linux (Ubuntu server 18.04), and, if possible, without passing through a PostGIS database?
A docker based solution could also be useful.
The tricky part will certainly be to set up a symbology (for which I don't know from where to explore at the time of writing...) but I wonder if one can use some predefined one such as the cartodb positron for example, or online tools such as https://cloud.maptiler.com/customize/#streets/0.82/0/0 or https://maputnik.github.io/editor/#0.83/0/0 to generate a custom style.json file which could then be used somewhere in a config?
What I've tried so far:
I've thought writing data to a sqlite3 database as follows could be a nice idea, but this will obviously not make it 'mbtiles' ready (i.e. it will not be encoded as protobuf chunks):
wget -O - http://www.example.com/some.pbf | ogr2ogr -f SQLite my.sqlite /vsistdin/
(this latter command is explained here)
I also gave GDAL's MVT driver a try, but this doesn't seem to work very well (I'm unable to serve the so build tiles, e.g. tileserver-gl complains that the tiles are not in the "openmaptiles" format (why?). And I cannot figure out if this format is implemented in GDAL. Maybe I need another server then?):
ogr2ogr -f MVT mytileset my_osm_data.osm.pbf -o METADATA_FILE='./tileset/metadata.json'
Creating MBtiles file from ".osm.pbf" file using TileMill is outdated for my usage.


texts: QGIS Put a label only on selected polygons

texts: I'd like to labelling only some polygons of a layer which are inside another polygon from another layer. Have you an idea how to do this? May be with a label rule?


texts: How to use field calculator to return the percentage of each record within a group (defined by a different field)? (ArcGIS Pro)

texts:
Within ArcGIS Pro, I am trying to summarize records in a shapefile by calculating a new attribute field that shows the percentage of a given record within each group, defined by another field in the attribute table.

Here is a dummy attribute data example of what I am trying to do. I am looking to use field calculator to calculate a new field GROUP_PCT:

ID, GROUP, VALUE, GROUP_PCT

1, A, 3, 75

2, A, 1, 25

3, B, 7, 43.75

4, B, 9, 56.25

So for record ID 3 (in group B):

GROUP_PCT = 7/(7+9)*100 = 43.75%

....record ID 3 makes up 43.75% of the sum of values in GROUP B.

Is there a smart simple way to do automate this with Python code block in field calculator (or otherwise)?


texts:
I'm using the OpenStreetMap plugin for QGIS and I've successfully rendered a .osm file (lines). I've tried a few styles from github but so far none of them label all of the streets like the osm default style.

Is the default style available somewhere? And what format should I be looking for for a .osm file, xml, sld, qml.


texts:
Working with a greycale scanned image: when I georef that image it rotates (expected) the blank corners are then filled with black. When the new image overlays I need to make those corners transparent. Obviously I cannot select black as a transparent Colour as the image is greyscale. Is it possible to persuade the GeoRef app to create those blank corners in, say, red - or some other colour so that I can make it transparent?

Allan


texts: Copying single selected polygon from large ArcSDE feature class in less than 10 seconds?

texts:
I have selected a single polygon, which is large but not massively complex, from an ArcSDE 9.3 for Oracle feature class with 3.5 million polygons using ArcGIS Desktop 10.0 SP3.

I would simply like to copy it out via ArcPy into a new shapefile (or file geodatabase feature class) and to do so in less than 10 seconds.

So far I have not been able to get it to happen in less than about 5-6 mins.

I think I need to try and take advantage of the processing environment (arcpy.env.extent) during the CopyFeatures, but if I set that to the values returned by getSelectedExtent and do the copy the time taken still seems to be 5-6 mins.

If instead I could get at a centroid x,y (and one that is guaranteed to fall inside the polygon) then I could use that to buffer a tiny amount and use that to set arcpy.env.extent and I would expect the time to come down to around 10 seconds.

However, to get at the centroid I think I still need to use CopyFeatures to get it out into a geometry first which seems Catch 22.

The 10 second target comes from Quickest way to select ArcSDE polygon by point in ArcGIS Desktop using ArcPy? which is another question based on the same dataset.


texts: Survey123 Connect geolocate from a text field using online survey

texts:
From the web survey, is it possible to put a point on the map from typing in your address? As it stands, when I type my address a drop down from the search bar in the map appears. But, the map doesn't zoom to my address.

Here's is my XLSX form:

From what I understand from this blog, I'm supposed to pull coordinates from the geopoint field, but I'm not sure if I'm doing this right. pulldata("@geopoint",${locate},"reversegeocode.address.Match_addr")


texts:
I've made some SVG symbols in Illustrator. Using them in QGIS 2.18 will give the elements in the SVG with no line, only a fill, on export to PDF automatically a line and fill.

SVG elements that have type in it, type also get on export to pdf a line and fill. Is this a bug or is there something I'm doing wrong?


texts: Using correct coordinate system in ArcMap?

texts:
Im doing fyp on sabah & saraawak map. it came without prj file. i manage projected it and its allign to each other. but now i think the map is too small. as in the whole continent is within few hundred meters.

What is wrong here?


texts: Layers not appearing even after source is loaded

texts:
I have written this code to source USPS postal route data from ESRI, which is to be displayed on a map.

The code works fine, the data source is also loaded, as per the console.log I put in place,  but I am not able to get the layers to be displayed.

Did I make a mistake in the layers block or otherwise?

The code is as follows.

let zipCode = 10001;
let esriAPI_1 = "https://gis.usps.com/arcgis/rest/services/EDDM/selectZIP/GPServer/routes/execute?f=json&env%3AoutSR=4326&ZIP=";
let esriAPI_2 = "&Rte_Box=R&UserName=EDDM";
let esriAPI = esriAPI_1.concat(zipCode, esriAPI_2);

var esrixhr = new XMLHttpRequest();
    esrixhr.withCredentials = false;
    //console.log("inside add event listnersasdads");
    esrixhr.open("GET", esriAPI);
    esrixhr.addEventListener("readystatechange", function () {
        if (this.readyState === 4 && this.status == 200) {
            ParseEsriData();
        }
    });
    esrixhr.send();

    function ParseEsriData() {
        esriData = esrixhr.responseText; //adding esri data to var
        geodata = JSON.parse(esriData); //parsing the data
        features_data = geodata.results[0].value.features; //returns the number of routes
        var eddm_parse = {
            'type': 'geojson',
            'data': {}
        };

        featlength = 1;
        for (let i = 0; i <= featlength; i++) {

            attribs = features_data[i].attributes;
            ckid = attribs.CRID_ID;

            path_length = features_data[i].geometry.paths;
            var geojson = {
                "type": "Feature",
                "features": []
            };

            for (let j = 0; j < path_length.length; j++) {
                //console.log(paths.length);
                //paths.length will give the total number of
                //path arrays from esri
                var marker = {
                    'type': 'Feature',
                    'id': ckid,
                    'geometry': {
                        'type': 'LineString',
                        'coordinates': ''
                    },
                    'properties': attribs
                };
                geocordinates = path_length[j];
                marker.geometry.coordinates = geocordinates;
                geojson.features.push(marker);
                marker.length = 0;
            }

            eddm_parse.data[i] = geojson;

        }
        console.log(eddm_parse);
        // map.addSource(ckid, { 'type': 'geojson', 'data': geojson });
        map.addSource('route', eddm_parse);
        map.addLayer({
            'id': 'route_layer',
            'type': 'line',
            'source': 'route',
            'layout': {
                'line-join': 'round',
                'line-cap': 'round'
            },
            'paint': {
                'line-color': '#888',
                'line-width': 8,
            }
        });

        function callthis() {
            if (map.getSource('route') && map.isSourceLoaded('route')) {
                console.log('source loaded!');
                map.on('data', function () {
                    console.log('route_layre');
                    map.getSource('route');
                    map.getLayer('route_layer');

                });
            }
            else {
                console.log("this didnt work");
                //return;
            }
        }

        map.on('sourcedata', callthis);


texts: How to create a custom enterprise QGIS installer with plugins and symbols?

texts:
I am looking at installing QGIS at a company level.  I would like my installer to install specific plugins, and add our company map symology.  Are there any ways to do this?  I have seen post about the custom GUI installer, but I didn't feel it was the same as my request.

Thanks...


texts: Print values from Tables in QGIS Console

texts:
I want to know how to get the QGIS Console to print out actual values pertaining to a Layer's attribute table (e.g. field names, feature values, etc.).

So far using the "print" command all I'm able to get are object references.

The next 2 examples involve a shapefile of mine called 'Province', and are predicated by the following code:

>>>vl = QgsMapLayerRegistry.instance().mapLayersByName('Province')[0]
>>>iface.setActiveLayer(vl)

Example 1:

>>>print vl

Output:

>qgis._core.QgsVectorLayer object at 0x11f24a348<

A vector layer object reference is returned, but I want the actual name of the layer, which of course is 'Province'.

Example 2:

>>>fields = vl.pendingFields()
>>>fields.toList()

Output:

[<qgis._core.QgsField object at 0x11f24a770>, <qgis._core.QgsField object at 0x11f24aa68>, <qgis._core.QgsField object at 0x11f24ab00>, <qgis._core.QgsField object at 0x11f24ab98>, <qgis._core.QgsField object at 0x11f24ac30>, <qgis._core.QgsField object at 0x11f24acc8>

Granted, I used "toList()" rather than "print", but the result is the same: rather than the actual names of the fields I get an array of object references.

How do you tell Console to give you the actual values that you happen to be storing in a given variable, rather than the object refs?

Much obliged.


texts: PyQGIS - Regex to Fill New Attribute Value

texts:
I have a gpkg with a column called name. I am trying to use PyQGIS to find any of the values in this column that starts with r and then contain e or x. I would then like to update a new column I've created called Relevant_Pipework with either Relevant or Not Relevant based on whether the regex has matched. Below is the code I've written.
import re
from qgis.core import (
    QgsVectorLayer
)

path_to_gpkg = "D:/Pipework.gpkg"

Water_Pipes_Layer = path_to_gpkg + "|Water_Pipes"

vlayer = QgsVectorLayer(Water_Pipes_Layer, "Water_Pipes", "ogr")
if not vlayer.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(vlayer)

layer = iface.activeLayer()
prov = layer.dataProvider()
pattern = '[rex]'

features = layer.getFeatures()

layer.startEditing()

if layer.dataProvider().fieldNameIndex("Relevant_Pipework") == -1:
    layer.dataProvider().addAttributes([QgsField("Relevant_Pipework", QVariant.String)])
    layer.updateFields()

for feat in features:
    Regex_Search_Field = feature['name']
    Regex_String = re.compile(pattern)
    Regex_Match = Regex_String.search(Regex_Search_Field)
    if Regex_Match:
        layer.changeAttributeValue(feat.id(), 11, "Match")
    else:
        layer.changeAttributeValue(feat.id(), 11, "No Match")


texts: What class/types of smartphone chipsets support GPS, A-GPS, GLONASS AND have extended kalman filtering?

texts:
I've been testing several Android phones (all running Android 4+)

Sony Xperia T, Moto G, Nexus and HTC 1 which are all giving decent tracks (all of which use GPS/GLONASS)

However we are also using the new Mi-Fone (Ai350) running android 4+ and we are seeing lots of noise. Specs are hard to come by for this phone but it looks like basic GPS. This is a sub $80 android phone

Question 1: What low end to mid-tier phones use GPS-GLONASS for better tracks?
We are using Google's 'My Tracks' app to record the data'

Question 2: Do GPS chips have extended kalman filtering on them? Some experts say yes while others say no. I am seeking some clarification on this.


texts: Removing the mirror route on PostgreSQL for shortestpath using pgr_ksp algorithm?

texts:
I have created an underground route for few points and when I try to see the possible routes using pgr_ksp I end up with an extra route which is not present in the geom table.
In the image below, this is one of the possible routes. In this route I am getting a line which is imaginary as there is no underground route to the diagonal.

How do I remove that line?
Below is the query used:
SELECT
  ST_AsText(ST_MakeLine(the_geom)),
  sum(cost),
  id1,
  sum(agg_cost)
FROM pgr_KSP(
       'SELECT structure_route_id::integer as id ,
         source::integer as source,
         target::integer as target,
         project_id,
         st_length(st_transform(the_geom,4326)::geography)::double precision as cost,
         st_length(st_transform(the_geom,4326)::geography)::double precision as reverse_cost,
         st_length(st_transform(the_geom,4326)::geography)::double precision as agg_cost
       FROM icw.underground_route
       where project_id=''"+ProjectId+"''
         and source is not null
         and target is not null
       order by source',
       "+source+",
       "+target+",10,false) as di
JOIN icw.underground_route pt
 ON di.id3 = pt.structure_route_id
group by id1
ORDER BY id1 ASC


texts: gdalwarp '-tps' equilvalent in C# wrapper?

texts:
I am attempting to do the following from code:
gdal_translate.exe --config GDAL_DATA "D:\gdal\bin\gdal-data" -of GTiff -a_srs EPSG:4326 -gcp 0 0 -108.724530873811 40.1376029956585 ..<truncated>.. -gcp 1920 1080 -108.700701203799 40.0986057211287 z.png y.tif

gdalwarp.exe --config GDAL_DATA "D:\gdal\bin\gdal-data" -of GTiff -wo "SAMPLE_GRID=YES" -wo "SAMPLE_STEPS=100" -overwrite -tps -r cubic -s_srs EPSG:4326 -t_srs EPSG:4326 -srcnodata "0 0 0" -dstnodata "0 0 0" y.tif y_done.tif

Starting with a non-georeference png file (obtained from a camera mounted on an aircraft), I add externally calculated GCPs to create a GeoTIFF.  Using that GeoTIFF, I want to warp that image to display it on a 2d map.
Using the -tps option in the gdalwarp command has so-far provided the best alignment for the resulting image.
In C# code, I am struggling to achieve the same final result using AutoCreateWarpedVRT:
Dataset dswarp = Gdal.AutoCreateWarpedVRT(dsTiff, s_srs_wkt, t_srs_wkt, ResampleAlg.GRA_CubicSpline, 0.125);
ds_out = dvTiff.CreateCopy(path + ".wrp.tif", dswarp, 0, null, null, "");

This creates a warped GeoTIFF but appears similar to an -order 2 (per the VRT file) polynomial rather than tps.
Am I going about this the right way?


texts: How to create sector figure in geoserver?

texts:
I have Azimuth, radius, angle, latitude and longitude.
I can create circular sector using polygon but it needs a lot of points and it will be slow to draw. Is there nicer way in GeoServer to do so?


texts: Changing the color of every entity of GeoJSON layer in OpenLayers 6.4.3

texts:
I would like to change the color of every entity of my GeoJSON, like in the picture. For now I can just change all the features with the same color.
Here is my code:
const provinceGeoJSON = new ol.layer.VectorImage({
          source: new ol.source.Vector({
              url: './data/chefchaouen.geojson',
              format: new ol.format.GeoJSON(),
              attributions: [ "&copy; INSEE", "&copy; IGN" ],

          }),

          title: 'provinceGeoJSON',
          style: function (f) {
        return [
          style,
          new ol.style.Style({
            stroke: new ol.style.Stroke({ width: 1.5, color: [0,0,0]}),

          })
        ]
      }


texts: How do you split multiple polylines by a field name?

texts:
I used the point to line tool in order to create lines going from one point to another. Each line essentially represents the route taken by an individual. How would I be able to separate the lines by the individuals name or objectid?

Each row has an independent line associated with it and I would like to take the shapefile shown in the screenshot and split it into multiple shapefiles depending on the name.


texts: Querying ARCGIS to find feature having specified point

texts:
I am trying to query features by point, so that it returns a feature that CONTAINS the point given. With the screenshot below, you can see that inspite of giving all the fields correctly, the server is returning all the features and not the ones which actually contains the given point.
(I am going to use the parameters later while querying through javascript api.)

I have tried after selecting all the spatial relationships but same result. Same for geometry type. I am not able to understand what I am doing wrong. I know that my database uses {wkid:32637} but with spatial reference given this shouldn't be a problem.


texts: On QGIS2.4 how can i represent a sens of circulation?

texts:
I'm working on the creation of a map of the third arrondissement of Lyon (France).
I'm using free informations from the web site "SmartData" of ("Le Grand Lyon).
I'm trying to represent roads, under QGIS 2.4, with arrows for the representation of the circulation (Oneway).
In my colonn i have 'FT' if it's in the sens of the circulation and 'TF' if it's on the other sens. So i want that the arrow go on inverse sens.
How can i do this with Qgis ?!


texts:
I have some queries which clip rasters from polygon.

So how I can see my result query in postgis or How I can check my query true or false ?


texts: How to generate WKT geometries from true type font glyphs (python)

texts: Does anyone know of a library or tool ( preferrably Python)to convert true type font glyphs into polygonal shapes (for testing GIS algorithms)?


texts: Making multivalue input script run from Python script tool?

texts:
I have a working python script:

import arcpy

sdeconnection = arcpy.GetParameterAsText(0)
layerfiles = arcpy.GetParameterAsText(1).split(";")

for layerfile in layerfiles:
    print layerfile
    lyr = arcpy.mapping.Layer(layerfile)
    lyr.replaceDataSource(sdeconnection, "SDE_WORKSPACE", lyr.datasetName, False)
    lyr.save()
    del lyr

I can run it from the command line with multiple layer files as input (separated by ";").

However I cannot make it work when I the script to add to a toolbox. The script runs as expected if I set "multivalue" to false for the layerfiles parameter.

However when I set "multivalue" to true I get the error:

Start Time: Thu May 16 10:13:47 2013
Running script changeLyrWorkspace...
<type 'exceptions.ValueError'>: Object: CreateObject Layer invalid data source
Failed to execute (changeLyrWorkspace).
Failed at Thu May 16 10:13:47 2013 (Elapsed Time: 0,00 seconds)

I get this error regardless of whether I add one or more layers when I run the script through the toolbox.

EDIT

It turns out that the script works through the GUI if there are no spaces in file name of the layer file.

Are there any workarounds for this problem?


texts: Why would geocoding by street address have ~96% failure rate?

texts:
I attempted to use CartoDBs Geolocate By Street Address functionality.

I used 3 columns: street address, City (Austin), and State (TX). I manually set the country to USA.

Of the 30 rows, only two ended up with a coordinate within the given city - the other 28 where all over the US. Is this typical?

Here is a link to the data: http://pastebin.com/w4MQgZuF


texts:
I'm looking for a high resolution raster image for the province of Maputo in Mozambique. I need it to get the elevation and topography for the area of Moamba-RessanoGarcía.
Anyone knows any website where can I find and download this raster file?


texts: Error Message Orthorectifying ASTER Imagery in ArcMap 10 with DEM

texts:
I am trying to orthorectify ASTER imagery (VNIR, 15m resolution) in ArcMap using an ASTER DEM  that I have downloaded from the GDEM website. The area is in the Andes, there are very high elevation differences so I guess orthorectification is needed.

I was using the tool "Create Ortho Corrected Raster Dataset (Data Management)" and I tried with the DEM as input and also with a derived DHM. I always get an error message very quickly.

It's the same coordinate system, I have chosen a very short path (C:/temp), I have tried with file formats *.img and *.tiff. I have never orthorectified before. Any ideas?


texts: R - Dealing with missing data in SpatialPolygonDataFrames for Moran.test

texts:
I'm trying to use Moran.test on a SpatialPolygonDataFrame consisting of 7194 elements in R. I know that there is around 150 polygons with NA values.

First I generate a spatial weights matrix:

W_cont_msoa <- poly2nb(all.data, queen=T)
 W_cont_msoa_mat <- nb2listw(W_cont_msoa, style="W", zero.policy=TRUE)

Then when I carry out Moran.test the code below returns:

> moran.test(residuals, listw=W_cont_msoa_mat, zero.policy=T)
Error in moran.test(residuals, listw = W_cont_msoa_mat, zero.policy = T) :
objects of different length

I've also tried to this but it hasn't worked:

> moran.test(residuals, listw=W_cont_msoa_mat, zero.policy=T, na.rm=TRUE)
Error in moran.test(residuals, listw = W_cont_msoa_mat, zero.policy = T,  :
unused argument(s) (na.rm = TRUE)

I've tried to go back to the original data and remove missing data using complete cases but I received the following:

> all.data[complete.cases(all.data)]
Error in all.data[complete.cases(all.data)] :
error in evaluating the argument 'i' in selecting a method for function '[': Error in
complete.cases(all.data) : invalid 'type' (S4) of argument

The structure of the data is as follows:

`Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
..@ data       :'data.frame':   7194 obs. of  23 variables:
.. ..$ MSOA01CD   : Factor w/ 7194 levels "E02000001","E02000002",..: 1 2 3 4 5 6 7 8    9 10 ...
.. ..$ MSOA01NM   : Factor w/ 7194 levels "Adur 001","Adur 002",..: 1468 125 126 127     128 129 130 131 132 133 ...
 .. ..$ averagesp  : num [1:7194] 11.6 20.9 20.3 21.4 19.1 ...
 .. ..$ nga        : int [1:7194] 41 112 177 97 120 NA 142 161 122 113 ...
 .. ..$ mediansp   : num [1:7194] 11.5 13.1 12.6 13.5 10 ...
 .. ..$ mbps2      : int [1:7194] 125 74 104 48 83 NA 64 67 76 25 ...
 .. ..$ msoa.total : int [1:7194] 2372 114 191 97 126 NA 157 208 130 119 ...
 .. ..$ p.nga      : num [1:7194] 1.73 98.25 92.67 100 95.24 ...
 .. ..$ p.mbps2    : num [1:7194] 5.27 64.91 54.45 49.48 65.87 ...
 .. ..$ residents  : Factor w/ 4057 levels "10,003","10,004",..: 2222 1748 27 1261       3216 NA 3380 403 3093 3256 ...
  .. ..$ area       : Factor w/ 6908 levels "1,002.23","1,002.58",..: 3707 2869 2840      3304 839 NA 1899 1327 6811 748 ...
 .. ..$ density    : num [1:7194] 25.5 31.3 46.9 24.8 72.1 NA 50.6 81.5 87.4 76.8 ...
  .. ..$ micro.bc   : num [1:7194] 36 32 83 33 39 NA 74 60 44 32 ...
  .. ..$ small.bc   : num [1:7194] 49 14 64 3 14 NA 50 14 14 73 ...
  .. ..$ medium.bc  : num [1:7194] 850 0 15 0 0 NA 10 0 0 0 ...
  .. ..$ large.bc   : num [1:7194] 245 0 0 0 0 NA 0 0 0 0 ...
   .. ..$ primary    : num [1:7194] 60 0 0 0 0 NA 0 0 0 0 ...
  .. ..$ secondary  : num [1:7194] 135 0 20 0 0 NA 0 0 5 5 ...
  .. ..$ tertiary   : num [1:7194] 16600 95 335 100 145 NA 300 200 160 80 ...
  .. ..$ kibs       : num [1:7194] 11805 20 80 20 30 ...
  .. ..$ edu        : num [1:7194] 2205 2681 811 2658 116 ...
  .. ..$ p.non.white: num [1:7194] 21.4 35 45.4 19 33.7 NA 32.8 30.4 38.8 30.8 ...
 .. ..$ claimants  : num [1:7194] 117 245 316 254 295 ...
 ..@ polygons   :List of 7194

Can anyone suggest how I overcome this problem please?

Thanks,

Ed


texts: Osm2pgsql importing speed not improving when using a SSD?

texts:
I have built my own OSM server and I am trying now to import UK into my Postgis database using osm2pgsql. I have changed the data directory of the database (succesfully) and I have mounted a Samsung 850 EVO SSD on it. The issue is, I do not see any speed increase when using SSD compared to a normal import using HDD.
Osm2pgsql configuration is like this:

osm2pgsql -d gis --create --slim  -G --hstore --tag-transform-script openstreetmap-carto.lua -C 3500 --number-processes 1 -S carto.style UK.osm.pbf

(Working in a Ubuntu 16.04 virtual machine of 5GB RAM)

Nodes and ways are being imported at a speed of 80 k/s and 0.10 k/s in average in both cases, respectively. Any idea about what could be happening?

*I am 99% sure I am importing into the SSD device as when I import a smaller region its free size decreases.


texts: Reclassifying slope using ArcGIS raster calculator?

texts:
I have a slope raster in degrees and would like to reclassify it the layer with the raster calculator.

my criteria: slope between 15 & 20 receive a 1
slope greater than 20 receive a value of 2
all other values receive a zero

Con(("noury_slope" >= 15) & ("noury_slope" <= 20,1,0)),Con("noury_slope" > 20,2,0)

get this error

ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 5, in rcexec
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\sa\Functions.py", line 3806, in BooleanAnd
    in_raster_or_constant2)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\sa\Functions.py", line 3803, in Wrapper
    ["BooleanAnd", in_raster_or_constant1, in_raster_or_constant2])
TypeError: expected a raster or layer name

Failed to execute (RasterCalculator).
Failed at Mon Aug 21 10:42:24 2017 (Elapsed Time: 0.30 seconds)

I want to do this query in one step. whats wrong with my syntax


texts: Search & Display Function for Geodjango-Leaflet Web Map

texts:
I am trying to add a search function into my Geodjango-Leaflet-Openstreetmap web map. I have around 540000 polygons of property boundaries that were drawn in Arcmap; the polygons will be imported to the web map later. I want to be able to search for each individual polygons by typing in the polygon's ID and have its exact position, polygon's shape, and attributes (address, dimensions, coordinates, owner, etc) show up on the web map.
I have just copied the code from here "leaflet search bar" into my web map and it seems to only work for markers that are predefined in the index.HTML page.
Is this something that is doable in Geodjango or am I better off using something else like ArcGIS Web Appbuilder?
Here is my .html code
<!DOCTYPE html>

{% load static %}
{% load leaflet_tags %}

    {% leaflet_js %}
    {% leaflet_css %}
    <title>Home</title>

    <link rel="stylesheet" type="text/css" href="{% static 'routing/leaflet-routing-machine.css' %}">
    <link rel="stylesheet" type="text/css" href="{% static 'leaflet-groupedlayercontrol/leaflet.groupedlayercontrol.css' %}">

Welcome

{% leaflet_map "gis" callback="window.map_layers" %}


texts: Converting multiband .nc files using R to upload to Google Earth Engine

texts:
I am new to using .nc files and am struggling to get the output I need in GeoTIFF format to upload to GEE.
I have downloaded the 2015 zip file of fire data from the CCI archive here: http://data.ceda.ac.uk/neodc/esacci/fire/data/burned_area/AVHRR-LTDR/grid/v1.0/Compressed. And have written a loop to load in the .nc files and export them as GeoTIFFs so I can use them in GEE. However, when I look at them in GEE they do not have band names. Secondly, I am expecting 11 bands as there are 11 variables detailed here: https://catalogue.ceda.ac.uk/uuid/4f377defc2454db9b2a6d032abfd0cbd. However, there are 18 in GEE.
So I am wondering if I have been able to load all of the bands in properly? (I need to use all of them).
Finally, I am not sure that a time variable is included in the data. Each .nc file represents a month, and I need to combine them all into one GEE image collection which I can filter by date.
My basic code is shown below.
### convert .nc files to .tif
### load libraries ###
library(ncdf4)
library(raster)
library(rgdal)

setwd("filepath\\Archive Fire CCI grid\\2015\\")`

### list all files in the folder
ncfiles <-  list.files()
ncfileTest <- brick(ncfiles[1])
ncfileTest
##loop over every file in the folder, exporting tif with simplified name
for(i in ncfiles){
ncfile <- brick(i)
## file name doesnt agree with GEE so needs to be cleaned
cleanname <- gsub("-","_", i)
cleanname2 <- gsub("_AVHRR_LTDR_fv1.0.nc","",cleanname)
writeRaster(ncfile, file.path("Output file path\\ESA CCI Fire Archive",cleanname2), format="GTiff",overwrite=TRUE)
}


texts:
Possible Duplicate:
  Raster generalization - buffers in rasters, expand pixels?

I have a 1m LIDAR DEM which I'm manipulating in QGIS.

I'd like to smooth the DEM, such as by having each pixel become the average of its 3x3 neighbourhood (a low-pass filter?).

Is there a way to do this in QGIS?


texts: Downloading and importing labels for CanVec Named Features (Code 1580010) layer (Natural Resources Canada)?

texts:
I used the CanVec symbology toolkit (described at Canvec Symbology in ArcGIS Desktop?) to symbolize all the downloaded data but now I am wanting to downland and import labels or a procedure that I should follow for creating all necessary labels.

I have an add-on that I can use to import labels from another layer so ideally I would like to obtain a layer file that already has all the labels classes defined and properly symbolized.


texts: How to disable zooming with scroll in OpenLayers2

texts:
I am trying to disable the scroll zoom functionality in OpenLayers.
I have tried different things but nothing works.
First I tried this:

map.addControl(new OpenLayers.Control.Navigation({'disableZoomWheel':true}));

But nothing happened.
Then I found this here and I tried it but didn't work either. Actually in this case the page doesn't even load anymore.

Any ideas what to do?


texts: How to transform ERA5 vertical datum to WGS84

texts:
I am trying to plot a netcdf file from ECMWF ERA5 in a folium map but I am having troubles with the projections in the file.
ECMWF documentation states:
ECMWF data is referenced in the horizontal with respect to the WGS84 ellipse (which defines the major/minor axes) but in the vertical it is referenced to the Geoid (EGM96).
The netcdf file in fact seems to line up well in the horizontal axis but not on the vertical axis.
Is it possible to convert the vertical projection to match the horizontal projection in WGS84?
Here is my python code:
import pandas as pd
from streamlit_folium import folium_static
import matplotlib.pyplot as plt
import folium
import numpy as np
import xarray as xr
import pyproj
from pyproj.crs import CompoundCRS

cmpd_crs = CompoundCRS(name="WGS 84 + EGM96 height", components=["EPSG:4326", "EPSG:5773"])
print(cmpd_crs)

ds = xr.open_dataset(r'/home/giacomo/anaconda3/envs/ERA5plot/2021_daily_max.nc')

ds1 = ds.reduce(np.nansum, 'expver')
print('netcdf opened')
print(ds)

sliced = ds1.sel(time=slice(pd.to_datetime('2021-01-28'), pd.to_datetime('2021-02-02')))
print(sliced)

map_data = sliced['i10fg'].max('time')
print(map_data)

np_image = map_data.values

### projection transformation
map_data.rio.set_crs(cmpd_crs)

inProj = pyproj.Proj(cmpd_crs)
outProj = pyproj.Proj('epsg:4326')

x_1, y_1 = map_data['longitude'], map_data['latitude']
xy = list(zip(x_1, y_1))
x, y = pyproj.transform(inProj, outProj, xy[0], xy[1])

### set coloration
val_min = float(np.amin(np_image))
val_max = float(np.amax(np_image))

cm = plt.cm.get_cmap('viridis')

np_image_norm = (np_image - val_min) / (val_max - val_min)

colored_data = cm(np_image_norm)

### create folium map
map_loc = folium.Map(location=[0, 0], zoom_start=0, crs='EPSG3857')
folium.raster_layers.ImageOverlay(
                name="Mercator projection SW",
                image= colored_data, #np_image,
                bounds=[[float(y_1.min()), float(x_1.min())], [float(y_1.max()), float(x_1.max())]],
                opacity=0.5,
                mercator_project = True,
            ).add_to(map_loc)
folium_static(map_loc)


texts: Linear Regression Google Earth Engine on seasonal dataset

texts:
How can I correctly display my slope in Google Earth Engine, at the Map.addLayer function? I don’t quite understand how to set up the min and max values on scale or offset. I want to show increasing, decreasing and no trends.

var addTime = function(image) {
  // Scale milliseconds by a large constant to avoid very small slopes
  // in the linear regression output.
  return image.addBands(image.metadata('system:time_start').divide(1e18));
};

var Zimbabwe = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw').filterMetadata('Country', 'equals', 'Zimbabwe');

Map.addLayer(Zimbabwe);

// Load several years of CHIRPS ( 2000 to 2017)

var collection = 'UCSB-CHG/CHIRPS/DAILY'; //daily precipitaion
var band = ('precipitation');
//var resolution = 500;
var startDate = '2000';
var endDate = '2017';
var startMonthDay = '-11-01';
var endMonthDay = '-04-30';

var startMonthDay2 = '-01-01';
var endMonthDay2 = '-05-01';

// Choose whether or not to add layers to map view and generate charts
var addToMap = 'TRUE';
var printCharts = 'TRUE';

var collection = ee.ImageCollection(collection).map(addTime);

Map.centerObject(Zimbabwe, 6);

if(addToMap=='TRUE') {
  var empty = ee.Image().byte();
  var regionBoundary = empty.paint({
  featureCollection: Zimbabwe,
  width: 3
});
  Map.addLayer(regionBoundary, {palette: '252525'});
}

// create a blank image ..read abt add.Band
var imageblank = ee.Image();

var list = [];
for (var i = startDate; i <= endDate; i++) {
  var iMod = i - 1;
  var collectionDates = collection.filterDate(iMod+startMonthDay, i+endMonthDay);
  var collectionMean = collectionDates.sum().clip(Zimbabwe);

  var collectionDates2 = collection.filterDate(i+startMonthDay2, i+endMonthDay2);
  var collectionMean = collectionDates2.sum().clip(Zimbabwe);

//create blank image and use .addband image to display the image collection

imageblank = imageblank.addBands(collectionMean);

Map.addLayer(collectionMean.select(['precipitation']),{min:400,max:1295,palette:['white','green','blue']}, 'Prec_Mean-before');

}

  if(printCharts=='TRUE') {
    var timeSeries = ee.ImageCollection(collection).select(band).filterDate(iMod+startMonthDay,i+endMonthDay);
    var display = {
    title: collection+' '+band+' '+iMod+startMonthDay+' through '+i+endMonthDay+' '+Zimbabwe,
    fontSize: 12,
    hAxis: {title: 'Year'},
    vAxis: {title: band},
    trendlines: {0: {color: '#aa3319', visibleInLegend: true}},
    series: {0: {color: '#161616'}}};
  print(ui.Chart.image.series(timeSeries, Zimbabwe, ee.Reducer.mean(), 1444).setOptions(display));
 }

print('imageblank',imageblank);

var linearFit = ee.ImageCollection(imageblank).select(['system:time_start', 'precipitation']).reduce(ee.Reducer.linearFit());

//Display the trend with increasing slopes in green, decreasing in red [bands:..].

Map.addLayer(linearFit,
  {min: 0, max: [-1, 8e-5, 1], bands: ['scale', 'offset', 'scale']}, 'fit');


texts: Spatial autocorrelation analysis (Global Moran`s I) suspicious z-score

texts: I have a doubt concerning the Spatial autocorrelation analysis Global Moran's I, in ArcGIS 10.2. I did the analysis with 3.372 features of a defined study area, and it showed that the features are clustered. But I am not convinced with the results, because they are too high: the p-value is 0, and the z-score, 40,12. Is this possible? The result seems suspicious.


texts: Saving Output of Algorithm using PyQGIS?

texts:
I'm new to PyQGIS and trying to save the output of the processing tool 'create points layer from table'. I just need to save it as a shapefile. My code is below which seems to be functioning so far. Using QGIS v3.4.

params = { 'INPUT' : 'C:/Users/xxx/WeeklyContacts-2019-07-15.csv', 'MFIELD' : None, 'OUTPUT' : 'memory:', 'TARGET_CRS' : QgsCoordinateReferenceSystem('EPSG:27700'), 'XFIELD' : 'Easting', 'YFIELD' : 'Northing', 'ZFIELD' : None }
processing.run("qgis:createpointslayerfromtable", params)


texts: Appropriate projection for US distances?

texts: I am generating distances between crash sites (points) and urban areas (shapes) and am wondering what the best projection is to use for this task.  The area of interest is limited to the U.S. (Alaska and Hawaii) included.  I could do Alaska and Hawaii separately to be more accurate if need be.


texts:
Is there a way to connect an excel table to ArcMap software or a geodatabase so that it updates when data in ArcMap is updated?
Or perhaps there is another software besides excel that could be used to create reports that will update.


texts: How do I convert WSG-84 coordinates to lat/long?

texts:
Possible Duplicate:
  coordinate convertion: decimal degree to deg-min-sec of  lat/long coordinate

I am looking for either an online tool or one that I can download to convert WSG-84 coordinates to lat/long that I am used to seeing.

For example, I want to convert from:

41.590833, -93.620833 to 41° 35′ 27″ N, 93° 37′ 15″ W ...  Those are the coordinates, according to a Geohack, for Des Moines.

Thank you.


texts:
I am a beginner in QGIS

I am planning to add line events on a route at 200m intervals.

I have done georeferencing of the map and digitized file in Kml format.

Can i convert this file into a shapefile?

I collected the road details of a 15 km road stretch at 200m intervals.

I want to display and analyse the collected features of the road in  GIS.


texts: Is it possible to replaceDataSource of layer with dataset stored in_memory using ArcPy?

texts:
I use in_memory workspaces to speed up some Python script tools that get published as Geoprocessing Services significantly.

The tools open MXDs and perform, among other things, an ExportToPDF.

In the MXD there is a layer that references a dataset stored on disk.

At runtime a new dataset gets created with the same schema as that above, so that it can replace it as the data source but, given a choice, I would like that new dataset to ONLY be written to the in_memory workspace.

My testing suggests that using in_memory workspaces like this is not supported by either Layer.replaceDataSource or Layer.findAndReplaceWorkspacePath.

Is there a way to replace the data source of a layer with a feature class stored in the in_memory workspace so that it can be displayed in a map that gets exported to PDF before the in_memory workspace is deleted?

It looks like this would possible using ArcObjects and an InMemoryWorkspaceFactory but here I need an ArcPy solution.

If there is no way to do this using ArcPy at ArcGIS 10.0 (and nothing coming at 10.1) then what I will do is to propose as an ArcGIS Idea that Layer.replaceDataSource be enhanced to support a workspace_type=INMEMORY_WORKSPACE option.

Note: I have also posted this to the ArcGIS Discussion Forums but have had no replies there so far.

An ArcGIS Idea has also been submitted for this to be enhanced.