texts: Create circles around a point on a map measuring the radius in hours

texts:
Disclaimer: I am not experienced in GIS and hope this question fits in this category.

I was wondering whether it is somehow possible to create a circle* around a city, using Google Maps or the like, whilst the circle is not measured in miles (or other distance measure) but in hours of driving.

The background is the following: I am looking for an ideal location to live if I have to commute between two different places. Ideally one place should be within a ~15 minute range and the other place in a ~45 minute range. I can't find an intersection manually. If there is an automated way how to do this, I could experiment with my parameters to find the ideal place.

Any ideas?

*Obviously it won't result in a circle.


texts:
I have a geotiff image that rings in at about 290MB. However, once I smooth it using a gaussian filter I'm unable to get it back to that size, it ends up being 1720MB. The size isn't really the problem so to speak, it's that when I go to use gdal_contour it takes 200x as long to contour the smoothed image.

The reason that I'm smoothing the image to begin with is that gdal_contour on the original DEM/DSM gives me a lot of little islands from noise in the data and very circuitous contour lines. Once I smooth it, the islands disappear and the contour lines are less wild, but the file gets so large that now the contouring takes forever.
here are my gdal commands:

Generate contours labeled with elevation at 1m intervals (takes about 90s)

    gdal_contour -a elev -i 1 original-dsm.tif

Generate contours labeled with elevation at 1m intervals (takes about 18000s)

    gdal_contour -a elev -i 1 smoothed-dsm.tif

I thought that maybe I could use

    gdal_translate (with various options such as -co COMPRESS=DEFLATE -ot Int16 and others)

to convert the smoothed tif to a smaller size and thus make the contouring go faster, but I've been unable to find the right combination of compression and bit depth.

I've used imagemagick to inspect the files and the major discrepancy I notice is that the original-dsm.tif file has "Channel depth: gray: 1-bit" whereas the smooth.tif has  "Channel depth: gray: 16-bit". I'm pretty sure that this is why it's taking so much longer to process and why the image file is so large. Another way I notice the manifestation of the same problem is that the Base Type of the original file (imagemagick identify command) is Bilevel, but the Base type of the smoothed tif is Grayscale.

Does anyone know how to change the channel depth using gdal_translate or something else? Keep in mind that this is georeferenced, and I can't afford to lose the georeferencing.


texts: Thiessen Polygon Total Area Calculations Model Builder/Python Code

texts:
I am attempting to create a tool using Model Builder (and Python code) for calculating the average annual precipitation of a study area using the Thiessen Polygon Method. My model runs good other than finding the total area of the study area for further calculations in my model.  Here is a sample of my code extracting the total area.  I attempting to just use Summary Statistics to find the total area but I don't know how to access that table for input into further calculations.

def MySub (TotalArea1) :

 TotalArea = 0.0
 cursor = arcpy.SearchCursor(feat)
 for row in cursor:
    TotalArea += row.getValue("POLY_AREA")
 del row
 del cursor

 cursor = arcpy.UpdateCursor(feat)
 for row in cursor:
    row.setValue("TotArea", TotalArea)
    cursor.updateRow(row)
 del row
 del cursor


texts: Trying to create a map with co-ordinates

texts: I'm trying to create a map with the co-ordinates for the location but if i select map units it creates black boxes in the corners of the map instead of a line of co-ordinates. If I used centimeters it puts the co-ordinates in the correct places but they are the same for latitude and longitude which is incorrect. I have ensured that the CRS is set correctly for the grid and the map and I'm struggling to work out how to rectify this. Below are screenshots to show the output I am getting.


texts: Raster to vector with gaps in Landsat ETM

texts:
I need to find size of water area from a raster (jpg) using QGIS 1.8.0-Lisboa. The raster is a Landsat image with gaps of data in it.

I used Polygonize, but the output file shows wrong polygon size exactly where the gaps are.

Is it because of the gaps in the image? Is there any other way to find size of water from the raster?

Also, after installing QGIS 1.8, Grass 6.4.3 doesn't open. Can not find an answer why?


texts: Highlight center polygon using Style OR Add Marker

texts:
We are working on Open source suit (Geo-server 2.8, PostgreSQL 9.6/PostGIS & Open Layer 3 etc.)

We developed web application using open layers 3 ,in which we are providing search tool (use to search plot number in map).

For ex. if we search '901' number then plot is coming on center part of the map.

Now we want to highlight OR place an marker on that for identification so we are checking this on our side and stuck at particular point.

We searched on Google related to this still not able to get an exact idea or some clue

Check this screenshot

following links/examples refereed on our side

https://jsfiddle.net/jonataswalker/sum62szv/

http://jsfiddle.net/expedio/6jtf8m5o/

http://jsfiddle.net/p_tsagkis/c4o4put8/

https://openlayersbook.github.io/ch10-openlayers-goes-mobile/example-02.html

https://openlayers.org/en/latest/examples/polygon-styles.html?q=custom

function plotZoomToExtend(geometry){
      //convert bounding box to points and calculate extent //
      var regExp = /\(([^)]+)\)/;
        var matches = regExp.exec(geometry);
        var seperator=" ";

        //matches[1] contains the value between the parentheses
        //console.log(matches);
        var res = matches[1].split(",");
        //console.log(res[0]);
        //console.log(res[1]);

        min_b=res[0].split(seperator);
        max_b=res[1].split(seperator);

        var extent=[min_b[0],min_b[1],max_b[0],max_b[1]];

        //console.log(extent);

        var extentStyle = new ol.style.Style(
        {
            stroke: new ol.style.Stroke({
            color: '#ffcc33',
            width: 3
        }),
            fill: new ol.style.Fill({
            color: 'rgba(0, 0, 255, 0.1)'
        }),
            text: new ol.style.Text({
            font: '12px Calibri,sans-serif',
            fill: new ol.style.Fill({color: '#000'}),
            stroke: new ol.style.Stroke({
            color: '#fff', width: 2
            }),
            text: ''
            })
        });

    // Create a polygon based on the array of coordinates
    var extentSource = new ol.source.Vector({});
    var polygon = new ol.geom.Polygon([extent]);

    // Add the polygon to the layer and style it
    var feature = new ol.Feature(polygon);
    feature.setStyle(extentStyle);
    extentSource.addFeature(feature);

    map.getView().fit(extent, map.getSize());
  }


texts: Reading Features From A Google Fusion Tables Table - OpenLayers

texts:
Why doesn't this example work? http://dev.openlayers.org/releases/OpenLayers-2.12/examples/fusiontables.html
I would like to do something similar but i can not find any other example how to display points from Google Fusion Table on the OpenLayers map?


texts: QGIS plugin requires a missing module: ui_gresp

texts:
I am developing my own plugin, the plugin manager shows the error:

  This plugin requires a missing module: ui_gresp

If i try to enable it the trace shows:

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

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 204, in     startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:\Users\rey/.qgis2/python/plugins\gresp\__init__.py", line 52, in classFactory
    from gresp import gresp
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\Users\rey/.qgis2/python/plugins\gresp\gresp.py", line 29, in
    from grespdialog import grespDialog
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\Users\rey/.qgis2/python/plugins\gresp\grespdialog.py", line 24, in
    from ui_gresp import Ui_gresp
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named ui_gresp

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.0.1-Dufour Dufour, d94c044

Python path: ['C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\\processing',         'C:/PROGRA~1/QGISDU~1/apps/qgis/./python', 'C:\\Users\\rey/.qgis2/python', 'C:\\Users\\rey/.qgis2/python/plugins', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGISDU~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISDU~1\\bin', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISDU~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISDU~1\\apps\\qgis\\python\\plugins\\fTools\\tools']

Ii'm not sure what this error means. I think it may be some error compiling the GUI.


texts: Changing pixel values of raster data from QGIS python plugin

texts:
I have been trying my hand at a QGIS plugin. The plug in reads raster layer data into an array, pipes the data to a C++ program, does processing over there, and pipes my results back. After a couple of weeks of fiddling with it, all that works fine.

I would now like to take these results and use them to change the pixel values on the processed raster layers. I have found numerous questions on the site that tell me how to read the raster data from the data provider, but I have not found one that explains to me how I can change the pixel values for the raster layer using Python w/ QGIS.

As a dummy example lets say I have an image of a forest and I want to change all the pixels in a band raster layer with values > 100 to -1 (for example).

How would I do that?


texts: How does grass r.viewshed take into account the Earth curvature?

texts:
Up to now I always had to do such analysis in a local Cartesian Coordinate System and on relatively small datasets, but today, I want to compute a viewshed analysis on a large raster DEM which is given in ESPG:4326 (WGS84) (lon, lat in degrees units) using GRASS GIS r.viewshed.
In the documentation of this function, I can read this for the flag -c:

-c
Consider the curvature of the earth (current ellipsoid)

I don't really now from that definition, if what they call the "current ellipsoid" is the one which is defined by the embedded projection of my loaded raster or the one defined when creating the grass project workspace (I'm totally new to that software)? This doesn't have much importance in my situation because both are the same but I'd like to know.
Then, I have the strong feeling that I may get into troubles, especially with long distances, with the fact that the DEM regular grid itself is given with (lon, lat) pixels expressed in degrees.
Hence my question;
Should I consider re-projecting the raster to a local cartographic plane or not, for both cases (with and without the -c flag) before running the viewshed analysis?

Projection type
-c flag enable
-c flag disabled

(initial data) WGS84 Geographic Coordinate System
A
B

Local Projected, Cartesian Coordinate System
C
D

For which of the A, B, C or D cases does or doesn't it make sense to run a viewshed analysis?


texts: Using Refactor fields (Batch processing): error in "Fields mapping" column

texts:
I have several layers and want to delete some of their table attributes columns, like FLAG. Using Refactor fields is fine:

However, when I do it in batch, I get an error in the Fields mapping column which doesn't recognize FLAG:

Missing parameter value: Fields mapping (row 1)

I was hoping to use the tool in batch as I have several columns I want to have deleted in many of my layers.


texts: OSM, Displaying a road with a line on a map from selected nodes

texts:
Along with a friend, we have started an educational side-project,for fun not as a commercial or trial software, to test what we have learned these last few months. We decided we would do a simple map navigation program that would take an input from a user (starting place + desired location) and draw the route based on the Open Street Maps data we have extracted.
For our first attempt we assume the user is in a car and implemented in Java a k-shortest path to find the route between them. I also used a short algorithm to find the nearest node to a user's present and desired location, which while not perfect, works neatly for the level we desire to achieve. And it's working really well, the routes given over to all test maps are always valid and generally fast, a few times even optimal.
The results are returned to an ArrayList of Nodes(id,ref,long,lat).

What we now want is to be able to input the starting-node,finishing-node and all the routing nodes as a single road in a map view in a web page.
Our general idea is that the user will be able to click two points at the map , and a red line will be drawn along the nodes-ways our Java class has returned.

How do you propose we do this?

I have tried leaflet and js, have only managed to center the user's desired map and zoom level, as well as a central point, but the line drawing has proven problematic for more than 2 nodes, either failing to display entirely or being severly out of road or in wild shapes. This has lead me to think I am taking an entirely wrong approach.

Could anyone please give me a suggestion as to what resources I should be looking at?
Thank you.

As for the code, I have to work off my head now. I seem to be missing something right now, but the general idea was as follows for the prototype:

var map;
var ajaxRequest;
var plotlist;
var plotlayers=[];

function initmap() {
    map = new L.Map('map');
    var osmUrl='http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png';
    var osmAttrib='Map data © OpenStreetMap contributors';
    var osm = new L.TileLayer(osmUrl, {minZoom: 12, maxZoom: 15, attribution: osmAttrib});

    map.setView(new L.LatLng(37.929, 23.6975),14);
    map.addLayer(osm);

    var roadLayer = new OpenLayers.Layer.Vector("Roads");

    map.addControl(new OpenLayers.Control.DrawFeature(roadLayer, OpenLayers.Handler.Path));

   var nodes = new Array(
      new OpenLayers.Geometry.Point(37.92, 23.69),
      new OpenLayers.Geometry.Point(37.93, 23.70),
      new OpenLayers.Geometry.Point(37.92, 23.70)
);

var roadLine= new OpenLayers.Geometry.LineString(nodes);

var roadStyle = {
  strokeColor: '#0000ff',
  strokeOpacity: 0.5,
  strokeWidth: 5
};
var drawnLine = new OpenLayers.Feature.Vector(roadLine, null, roadStyle);
lineLayer.addFeatures([drawnLine]);
map.addLayer(roadLayer);
}

The array of Geometry points is the Array that is returned by the Algorithm, but for the example I used a few static points. In the actual code, the line is returned severly out of bounds for a few nodes, while it is not returned at all after the nodes increase. However for 2-3 nodes it is working correctly. That's why I assume my approach is bad.


texts: Converting Lon Lat to UTM Zone 32N in C# with SQLSpatialTools

texts:
I have a set of coordinates; Lon 9.42221521, Lat 56.44978001 and I want to convert them to UTM Zone 32N coordinates. I found the SQLSpatialTools library to do the conversion but I need some help. When I use the following code:

SqlChars chars = new SqlChars(string.Format("POINT({0} {1})", lon, lat));
        SqlGeography mapPoint = SqlGeography.STPointFromText(chars, 4326);

        SqlProjection proj = new SqlProjection();
        proj = SqlProjection.TranverseMercator(0);  //(or Mercator)
        SqlGeometry geometry = proj.Project(mapPoint);
        geometry.STSrid = 25832;

        SQLSpatialTools.AffineTransform affTrans = new AffineTransform();
        SqlGeometry newPoint = AffineTransform.Scale(6378137, 6378137).Apply(geometry);

...the lon lat set gets converted to 578664.03 - 6323756.73

If I check the lon lat set in ArcMAP I get the following values: 526026.284 - 6256221.368

...so obviously I must be doing something wrong. I'm using SQLSpatialTools because it's easier for me to perform some spatial queries later.

Thanks!


texts: Is there any way to use mbutil on Windows?

texts: I have used mbutil on Ubuntu. It's a very handy tool. I am curious if there is any way to use it on a Windows 64 bit machine. If there is a way please help me to do so.


texts:
I've got the following task to solve:

I've got two tables; one contains buildings (houses) as points and the other contains train-stations (at_bhf_classified) as points.
I need to find out the distance from each house to the next trainstation.
Well, this is more or less easy.
But now the issue: The trainstation are categorized. That means:
Give the distance to the next station of category 1 if within a 25km radius. If no category 1 station is found within this distance, give me the distance to cat 2 station within the 25 km radius.

So far I use the following statement to calculate the nearest category 1 station:

WITH dist_min
 AS (SELECT h.h_key,
            Min(St_distance_sphere(h.geom, b.geom)) dist_min
     FROM   house h,
            at_bhf_classified b
     WHERE  St_distance_sphere(h.geom, b.geom) < 25000 and b.kat=1
     GROUP  BY h.h_key)
SELECT h.h_key,
   b.bhfname,
   b.kat,
   b.bhfnummer,
   round(cast(d.dist_min as numeric),2)
FROM   dist_min d,
   house h,
   at_bhf_classified b
WHERE  St_distance_sphere(h.geom, b.geom) = dist_min
   AND St_distance_sphere(h.geom, b.geom) < 25000
GROUP BY h.h_key, b.bhfname, b.kat, b.bhfnummer, d.dist_min
ORDER  BY h.h_key

Do you know how to include the second condition explained further up?

Thanks for you help!
Carsten


texts: Dividing polygon and assigning values proportionally (population) using QGIS?

texts:
I need to split several polygons (census blocks)  that contains numerical attributes (population) into smaller blocks, but assigning that data proportionally.

For example: if the splitting is made by 40/60%, and the population value is 1000, i want QGIS 3.2 to assign 400/600 proportionally to each polygon.

With the Split features function, the data is duplicated (QGIS replicates ALL the original values). I would love if there was a plugin that allows it to be made using the Split Polygon function "by hand", but I know that it is practically impossible.

The calculated field does not work.

1st step: create a calculated field ($area)

2nd step: cut the “one” polygon:

3rd step: create a calculated field (“population” * $area / “old_area”) as Virtual Field:

Result: the first result is OK (on the the polygon that was cut first), the “new_pop” value was calculated fine:

BUT the next split does not work (polygon "two"), the “new_pop” values are: 100 and 40, when should be 60 and 40. The first row is not updated.

Is this a bug?


texts: Is it possible to edit QGIS features in a web client by using QGIS server?

texts: I am new to QGIS and I am trying to build an application which runs on the web by using QGIS server. Is it possible to edit features or add features on web by using QGIS server?


texts:
So I'm converting a set of polygons into rasters and need to remove a portion of their basenames while keeping them unique at the same time. The error code I am getting is "TypeError: unsupported operand type(s) for -: 'unicode' and 'str'" obviously it has something to do with the " - ". I have tried to designate "_buffer" as a string [str("_buffer")] and that didn't work either.

here is the code

import arcpy
from arcpy import env
env.workspace = r"N:\geog327\python15\Project"
env.overwriteOutput = True

DenMetHwyBuff = r"N:\geog327\python15\Project\DenMetHwys_buffer.shp"
DenMetIntBuff = r"N:\geog327\python15\Project\DenMetInts_buffer.shp"
DenMetRiverBuff = r"N:\geog327\python15\Project\DenMetRivers_buffer.shp"
DenMetLakeBuff = r"N:\geog327\python15\Project\DenMetLakes_buffer.shp"

vectorlist4 = [DenMetHwyBuff, DenMetIntBuff, DenMetRiverBuff,
DenMetLakeBuff]
outraster = desc.basename - "_buffer.shp" + ".tif"
for fc in vectorlist4:
    desc = arcpy.Describe(fc)
    outraster = desc.basename - "_buffer.shp" + ".tif"
    arcpy.PolygonToRaster_conversion (fc, "Bad", outraster)


texts: Planet filter file for ENVI atmospheric correction

texts:
I have 4 band Planet Scope multispectral imagery that I am trying to process in ENVI.  It is the analytic data, so it is already radiometrically corrected.  I have converted the .tif to a BIL, and done the band specific calculations to get TOA reflectance (from the metadata).  I want to do the FLAASH atmospheric correction, but in the hyperspectral tab of the FLAASH window, it requires a filter function file that is specific to the sensor being used. Usually, with a known sensor, this filter function file gets automatically uploaded because ENVI saves .hdr and .sli files of many known sensors. I have ENVI 5.3 (and dont have the capabilities to upgrade).  In ENVI 5.5 they include the Planet Scope sensor and its filter details.

Does anyone know how to create a filter file for planet scope, and/or the format required for the file?


texts: Uncaught SyntaxError: Unexpected token : when using ajax

texts:
I try to call json data from geoserver,
but when i try to run it in chrome, it return this following error

"Uncaught SyntaxError: Unexpected token :"

This is my code,

    function initialize(){var urljson='http://localhost:1234/geoserver/Ven/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Ven:alamatpuskesmas&maxFeatures=50&outputFormat=json&format_options=callback:getJson';

             $.ajax({
                    url :urljson,
                    dataType: 'jsonp',
                    jsonpCallback: 'getJson',
                    success: handleJson
                });

        }
        function handleJson(data)
        {
            console.log(data);
        }

i am using geoserver 2.4.2, whats wrong with my code?,
when i try to paste the url in chrome/firefox, it return json. .


texts:
This is my website for test: http://www.cv.nctu.edu.tw/shapefile-js-openlayers/ol_simple.html

My dem data is (.jpg).

var inProj = new OpenLayers.Projection('EPSG:3826');
var outProj = map.getProjectionObject();
var options = {
    opacity: 1.0,
    isBaseLayer: false,
    numZoomLevels: 20
};

var extent = new OpenLayers.Bounds(16120,2409591,484824,2814500).transform(inProj,outProj);

var imageLayer = new OpenLayers.Layer.Image('Image Layer',
    'http://www.cv.nctu.edu.tw/shapefile-js-openlayers/naturalearthdata/cultural/dem.jpg',
    extent,
    new OpenLayers.Size(1,1),
    options);
map.addLayers([osm, GridshpLayer,BounderyshpLayer,RivershpLayer,SeashpLayer,imageLayer]);

error:


texts:
In ArcGIS Desktop 10.7 we have IDW tool, but when you want to work with it you should need two layers such as : Some points and a closed polyline (Boundary) as following:

then using IDW:

but I have a problem. my data is not point is a polyline (Actually is a line joining points of equal rainfall on a map in a given period) as following:

How should I use IDW tool?

This is because I don't have the lines in the first box of IDW tool:


texts: Do all open source web mapping platforms require programming?

texts:
I currently use a proprietary software product for my company's web mapping, however I am interested in seeing other products. I use QGIS and would like to to put my layers on the web for both internal (very large files - all roads in UK; all rivers in UK) and client based (small - 5-10 feature point/polygon files). From what I gather I need to use another platform (Leaflet, OpenLayes) in order to combine my layers into an online map with viewable attribute tables/searchable postcodes.

Do all of these platforms require JavaScript knowledge?  It seems like they all require coding of some sort and I cannot just 'upload' my shp files to an online server. If this is true, is one 'easier' than the rest?


texts: Error in raster displays of NoData cells when "zoomed in" past scale 1:3,000,000?

texts:
I have encountered an error in how ArcGIS Desktop displays a Raster file (.TIF format) at different scales. At scales greater than 1:3,76,600, when I'm zoomed in, the raster displays mostly correctly, showing only a border of low-value cells, or no border -the "border" appears inconsistently. (first picture)

What is consistent is that when I zoom out, specifically when I go from 1: 3,000,000 to 1:4,000,000, the white background (correctly indicating no-data cells) change to blue (indicating the lowest value "bin" in the classification). When I check the cell values, the cells are NoData value, even when zoomed in and blue (and they should be clear).  (second picture)

Any suggestions?

I want to display this data zoomed out to a scale smaller than 1:3,000,000 but it is falsely visualizing the data. I am working in ArcGIS Desktop 10.2 with all licenses installed. In Layer Properties > General Tab > "Scale Range," "Show layer at all scales" is checked.


texts: Is there R method similar to ArcGIS StandardDistance_stats?

texts:
I am finding some methods that can provide simple output.

I found a useful method StandardDistance_stats in ARCGIS, but the result is not what i am looking for, I only need the result of a lat/lng point and the size of the circle in number form, so I can plot it to my web map.

Now I am using geomean in geosphere package, it gives the mean of lat/lng points but still missing the size of the circle.

Edit: I found a solution for it, but I am stuck in import R script into Java

d <- read.csv(file="./points.csv", header=TRUE, sep=",")
x <- d$lng
y <- d$lat
w <- d$weight
df <- data.frame(x,y)
mean = geomean(df,w)
.sdd <- function(t,mean) {
    x <- t$lng
    y <- t$lat
    x <- x-mean[1]
    y <- y-mean[2]
    x <- x*x
    y <- y*y
    result <- (Reduce("+",x)/length(x))+ (Reduce("+",y)/length(y))
    result <- sqrt(result)
    return(result)
}
.sdd(d,mean)

Edit: I was trying to run R script file in java. i figured it out. It was because the permission problem, I was not allowed to access the root of C drive.

Base on this answer, I am not sure what should I modify in this line. I know that there are lots of Datums. if my list of lat/lng are mostly 22.0,144.0m is this line still suitable?

CRS.new <- CRS("+proj=somerc +lat_0=46.9524056 +lon_0=7.43958333 +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs")


texts: Using WCS and Python to download from USDA Vegscape database

texts:
I wish to download Weekly NDVI rasters from the USDA Vegscape portal. It is possible to access the data using WMS and WCS methods. The USDA dev guide only contains example URLs on how to access the data, see below:

https://nassgeodata.gmu.edu/VegScapeService/weekly_ndvi_2012.cgi?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&coverage=Weekly_NDVI_30_2012.07.24_2012.07.30&CRS=epsg:102004&bbox=130783,2203171,153923,2217961&resx=250&resy=250&format=gtiff

I found the below section of code on gis.se but am unsure how to adapt it for my own purpose:

import os
from owslib.wcs import WebCoverageService as w  #owslib ->https://github.com/geopython/OWSLib
                                                #easy_install owslib

folder = 'c:\\path\\to\\folder\\'
wcs = w(url="http://dmcsee.org/cgi-bin/mapserv?map=/var/www/tmp/dmcsee_wms/dmcsee_wcs.map",version="1.0.0")
spi6 = w['SPI6']  #The layer I am interested in

for timep in spi6.timepositions:
   output = wcs.getCoverage('SPI6',time=[timep],bbox=(10.0,31.995,48.005,50),format='GTiff',CRS='EPSG:4326', WIDTH=380, HEIGHT=120)
   f = open(os.path.join(folder,'spi6_'+str(timep)+'.tiff'),'wb')
   f.write(output.read())
   f.close()


texts: Upstream trace in geometric network by number of segments

texts: I am using a geometric network of a sanitary sewer. I am trying to determine the source of grease in the individual sewer segments. I know which segments that have been impacted by grease. I am trying to find the next 10 upstream segments.  In a geometric network is there a way to conduct an upstream trace that only goes up 10 (or what ever number I choose) segments from the segment where the trace started?  I am not looking for a distance, I am looking to identify specific upstream segments. I am using ArcMap 10.2.2.


texts: Tile services ON TOP of default CartoDB Basemaps?

texts:
I'd like to place a tile service of a georeferenced map from DavidRumsey.com in a CartoDB map, but I don't want it to be the 'basemap'. I want it to lay on TOP of the basemap, which would be any of the default choices.

Any way to do this within the GUI? Obviously it can be done with Leaflet, etc... but I'm writing guidance for a non technical class to use CartoDB for a History of San Francisco project, so I'd like to keep them in the GUI.


texts: QGIS 3.4.4 and 3.6 SAGA Tools all not working with same error

texts:
C:\OSGeo4W\bin>exit
Execution completed in 12.17 seconds
Results:
{'SHADE': 'C:/QGIS/output/HS.sdat'}

Loading resulting layers
The following layers were not correctly generated.C:/QGIS/output/HS.sdatYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

I get the above message on all the SAGA tools that I'm trying to run in QGIS, could someone give me an idea of what to look for here? I've tried switching versions of the ltr in the OSGeo4W installer, same thing happens.

Update (from comments):

I don't need help with any specific dataset per se, I think something is wrong with the way I am installing (64-bit advanced mode network installer OSGeo4W).

I've tried several times to install the 3 pieces, I need each of them, to no avail. There is always errors, i.e. the path variables written, versions, location of GRASS "C:/OSGeo4W/bin/bin" (notice two bins there) gets written and then GRASS doesn't work.

I was only able to get GRASS and SAGA both working by installing the 64-bit standalone installer long term release. All combinations of OSGeo4W installer, various ltr, latest, etc., all wrote either bad PATH variables, or one or more tools were pointing to the wrong folder locations for binaries. Which might not be a problem except I can't seem to find the processing option/dialogue where I can tell QGIS where the GRASS folder is.

And this is after getting a third party uninstall tool to clean up Windows registry and deleting my user profiles, otherwise problems would come back.

I'm documenting this in case other Windows users are in the same boat.


texts: Dynamically change location marker based on user's input

texts:
I have the following code for displaying OSM with a predefined feature. I am using OpenLayers3:

var wmsLayer = new ol.layer.Tile({
            source: new ol.source.OSM()
        });

        var view = new ol.View({
            center: [637125.42195, 8172199.19090669],
            zoom: 14
        });
        var iconFeature = new ol.Feature({
            geometry: new ol.geom.Point([637125.42195, 8172199.19090669]),
            name: 'Null Island',
            population: 4000,
            rainfall: 500
        });

        var iconStyle = new ol.style.Style({
            image: new ol.style.Icon(({
                anchor: [0.5, 46],
                anchorXUnits: 'fraction',
                anchorYUnits: 'pixels',
                opacity: 0.75,
                src: 'https://evernote.com/media/img/getting_started/skitch/android/android-location_icon.png'
            }))
        });

        iconFeature.setStyle(iconStyle);

        var vectorSource = new ol.source.Vector({
            features: [iconFeature]
        });

        var vectorLayer = new ol.layer.Vector({
            source: vectorSource
        });
        var map = new ol.Map({
            layers: [wmsLayer, vectorLayer],
            target: 'map',
            view: view,
            //controls: []
        });
        map.on('singleclick', function (evt) {
            var position = evt.coordinate;
        });

How I can change the location of the marker in the click event now?


texts: Hex grid on EPSG:4326 globe (for use in Leaflet)

texts:
I'm looking for a way to evenly distribute a hexagon grid (x,y,z) on the EPSG:4326 globe. Knowing this is not mathematically possible I'm accepting the following limitations:

Poles don't need coverage - from  lat -75 to lat 75 is sufficient.
International date line - 0/180 lng is OK to be messed up.
Distortion is OK as long as its not noticed by eye when viewing the grid at zoom ~100*100 hexes.

My hexagonal grids topology is defined according to this picture:

(source: keekerdc.com)

Some further conditions

The hexagons are ~50 meter height.
I'm assuming hex(0,0,0) to be located at LatLng(0,0)

I need a function hexToLatLng(x,y,z) and a function LatLngToHex(lat,lng)
How can these lookup functions be created in JavaScript?

Below are some of my attempts.

By using the estimates

M meter in north direction = M *(1/111111) degrees lat
M meter in east direction = M *(1/(111111*cos(lat)) degrees lng

I implemented hexToLatLng and LatLngToHex. However increasing errors seems to sum up and give big shifts in longitude accuracy when moving away from equator. I guess its because the estimates assuming a sphere and EPSG:4326 is ellipse?

Because the hex system is planar I first tried EPSG:900913 and proj4js until I learned EPSG:900913 is not very suitable for distance calculations. I would also prefer sticking to EPSG:4326 as its easy to work with in my Leaflet map. But any thoughts on this track also most welcome.

I did find this thread but I'm hoping there are other solutions given I'm accepting the drawbacks at the poles
https://stackoverflow.com/questions/749264/covering-earth-with-hexagonal-map-tiles


texts: GDAL tool for converting raster to point shapefile based on locations of another shapefile?

texts:
There is a shapefile with point data of rainfall sampled at a regular grid, and there is a shapefile with point data of soil attributes sampled at a regular grid, but the latter with smaller pixel size (finer grid). For each point on the finer grid I need to obtain an interpolated value of rainfall data. Basically, I was thinking of converting rainfall data into a raster using gdal_grid and interpolating it using IDW method.

Are there any GDAL tools which will enable to convert the raster back to a shapefile with point data matching locations of the soil properties (given in shapefile)?


texts: Error in ArcGIS Pro when Adding Layer to Group through ArcPy

texts:
I'm trying to add layers to an empty group layer through arcpy in ArcGIS Pro. I've read that creating empty group layers only works by importing a saved template. I'm able to load my template and it appears in the contents pane. But when I try to add a layer I get the following error.

# Import template
dir_path = os.path.dirname(os.path.dirname(__file__))
in_path = os.path.join(dir_path, r'templates\group_template.lyrx')
groupLayer = arcpy.mp.LayerFile(in_path).listLayers()[0]

# Add group layer to map
aprx = arcpy.mp.ArcGISProject("CURRENT")
aprxMap = aprx.activeMap
aprxMap.addLayer(groupLayer, 'TOP')

# Get layer to add
lyr = [l for l in aprxMap.listLayers() if l.name == 'test'][0]

# Add layer to group
aprxMap.addLayerToGroup(groupLayer, lyr)

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\utils.py", line 191, in fn_
    return fn(*args, **kw)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 1317, in addLayerToGroup
    return convertArcObjectToPythonObject(self._arc_object.addLayerToGroup(*gp_fixargs((target_group_layer, add_layer_or_layerfile, add_position), True)))
ValueError: <MappingLayerObject object at 0x0000000044B73770>

If I add the layer with aprxMap.addLayer(lyr) everything works fine and groupLayer.isGroupLayer returns True (lyr is a feature layer). Also, I'm able to add the layer by drag and drop, but I don't want to do that, since this code should finally be part of a larger function. I also tried converting the layer to a layer file as proposed in this question Adding feature class to map (within specific group) using ArcPy with ArcGIS Pro?, but I still get the error.

Does anyone know what I'm doing wrong?


texts: Create DEM with multi resolution and multi scale

texts:
How can I create a DEM with multi resolution and multi scale
I want to change DEM of details when the scale changes
and also I need high resolution DEM


texts:
I have been trying to change the height of my leaflet map to a percentage inside of bootstrap but everytime I do the map will not draw. Thus, I always have to revert to px value. I'm pretty sure its a simple setting that I'm missing since I am a CSS novice. Here is my css.

<style type="text/css">
  body {
    padding-top: 60px;
    padding-bottom: 40px;
  }
  #map {
    height: 75%;
  }
</style>


texts: How to mirror two geodatabases on the fly?

texts:
I'm new to arcgis and its applications.

We have the full suite of arcgis solutions (arcgis server, desktop, arcgis sde)

We have a requirement...

geodatabase (A & B)
A & B reside on 2 different server but connected via LAN
A will be updated
B will be a replica of A

May I know what is the best solution to this issue? the geodatabase can be hosted on oracle DB, personal geodatabase or geodatabase file.

I know that arcgistools has the functions to replicate and sync but according to my friend. the synchronization has to been done manually. (creeate a script and execute it periodically.

Is there a way to mirror those 2 geodatabase on the fly?


texts: how to enable interactivity on named map in cartodb

texts:
I have created a named map that gets data from a private table.  I can load it into a Leaflet map and can change its marker-fill color placeholder using layer.setParams.

However, I cannot get interactivity to work.  I want to be able to hover or click on a feature and do something useful.[1]

What do I need to do to make hover or click events fire?

Here is the named map json file

{
  "version": "0.0.1",
  "name": "mynamedmap",
  "auth": {
   "method": "open"
  },
  "placeholders": {
      "color": {
        "type":"css_color",
        "default":"blue"
      }
  },
  "layergroup": {
    "version": "1.0.1",
    "layers": [{
      "type": "cartodb",
      "options": {
        "cartocss_version": "2.1.1",
        "cartocss": "#layer { marker-fill: <%= color %>; }",
        "sql": "select * from mytable",
        "interactivity": ["cartodb_id"]
      }
    }]
  }
}

Here is the javascript to create the layer:

cartodb.createLayer(map, {
        type: 'namedmap',
        user_name: 'myusername',
        options: {
          named_map: {
            name: 'myusername@mynamedmap'
          }
}})
.addTo(map)
.done(function(layer) {
    console.log(layer.getSubLayerCount()); //prints '0'
    layer.setParams('color', 'green'); // works to change marker color

    layer.on('featureClick', function(e, latlng, pos, data) {
      console.log(e, latlng, pos, data); //this never gets called
    });

  });

[1] "Something useful" could be as simple as open infowindow and display cartodb_id and feature name.  Or it could mean "get cartodb_id for feature, pass it to my proxy where I can use curl to run some custom sql against my private table, then return json to the map and update display"


texts: Migration plugin from QGIS 2 to QGIS 3: list layers

texts:
I am trying to migrate a plugin from QGIS 2 to QGIS 3. This function doesn't work with QGIS 3 :

def getLayersByKeyWord(keyword):
    #retourne une liste de layer
    res=[]
    layerMap = QgsMapLayerRegistry.instance().mapLayers()
    for name, layer in layerMap.iteritems():
        if layer.keywordList() == keyword:
            res.append(layer)
    if len(res)>0:
        return res
    else:
        return None

I tried to replace QgsMapLayerRegistry with QgsProject but it is still not working on iteritems().

I still have this error :
for name, layer in layerMap.items:
TypeError: 'builtin_function_or_method' object is not iterable


texts:
Drawing ellipse with shapely?
How do I draw an ellipse with shapely?


texts: Arcade-get 4 nearest addresses from point feature

texts:
I used the below code to get the nearest one location from $feature, but I'm looking now to get nearest 4 records from $feature.

I tried with json object, arrayList, but having syntax errors. Not sure about the limitations of customization in Arcade.

var searchDist = 6000;
var storeAddress = Intersects(FeatureSetById($map, "storeAddressFeaturelayer_2246"), Buffer( $feature, searchDist, "meters"));
var cnt = Count(storeAddress);
var x=0;
var nearestaddress;
var neardistance=0;
var addressDist=0;

for (var f in storeAddress {
  addressDist = Round(Distance(f, $feature, "meters"),3);
 if(x==0){
     neardistance=addressDist;
     nearestaddress = f.Address+",  "+f.PhoneNumber;
  }
 if(neardistance>addressDist && x==1){
     neardistance=addressDist;
     nearestaddress = f.Address +",  "+f.PhoneNumber;
  }
  x=1;
 }

return nearestaddress;


texts: Rendering small polygons on high zoom using GeoServer?

texts:
I've noticed that GeoServer 2.7.6 has the capability to NOT render on high zoom small polygons. And GeoServer 2.11 is rendering all the features/polygons from a layer.
How can I make GeoServer 2.11 to stop rendering small polygons on high zoom?!
I've searched answers on Google but I didn't find info about this topic.
Only SLD scale dependent rule.


texts: How do I find forest areas over 1,000,000 acres with no roads?

texts:
I'm using QGIS and would like to identify all areas where intact forest is greater than 1,000,000 acres and not impacted by roads. In other words I need to identify polygons 1m acres + with no roads inside.

I have two global dataset (both vectors) - one of forest cover and one of roads. I would like to have a set of polygons showing these intact forest patches.


texts: Selecting layer by default using Leaflet?

texts:
I am working with the leaflet, where we first load the imageoberlay to the map using user input.

Now, i been trying different methods to select that layer by default.

The problem is only that the first layer not loaded on map until i click on the radio button.

How can i select layer by default. I am doing something like :

Script

        // New layergroup, note it's not added to the map yet
        var layerGroup = new L.LayerGroup(),
            imageOverlayUrl = firstImage,
            // New imageoverlay added to the layergroup
            imageOverlay = new L.ImageOverlay(imageOverlayUrl, bounds).addTo(layerGroup),
            // New featuregroup added to the layergroup
            featureGroup = new L.FeatureGroup().addTo(layerGroup);

        //dynamic method for adding new layer
        var layerGroupings = { "Main": layerGroup };
        var layerControl = new L.control.layers(layerGroupings, null, { collapsed: false });
        layerControl.addTo(map);

When new layer loaded, it should be selected by default.

How can i do that?


texts:
I can't see whether Esri CityEngine supports VRML files directly, or how to load them if they're not supported (the online materials are more marketing-focussed than technical).

Is there a recommended workflow to bring VRML files into CityEngine?


texts: How to create NDVI map From LiDAR intensity and RGB image

texts:
I am new with GIS and I've been looking for ways on how I could create a NDVI map from LiDAR intensity and RGB image. Is there any way I could do this?

I have read that you could create an intensity image but then how can we merge it with the RGB image? Is there a GIS software for that?


texts: How to install sda4pp plugin with missing repository in QGIS

texts:
I see that the repository for the sda4pp plugin went down in 2011.  Does anyone know if there is another way to download the plug-in with the repository missing?

If this plug-in is no longer available, does anyone have recommendations for kriging?


texts: Getting raster properties Using Mask File in ArcPy?

texts:
I am trying to get raster properties using a Mask file for pixels that fall within mask, script will calculate that area statistics only.

Below is my script I am working on, but its showing an error

import arcpy
arcpy.env.workspace = r"D:\Annomaly_SPEI_12"
Mask = r"D:\Agro_Zone\Zone-1.shp"
rasterList = arcpy.ListRasters()

for raster in rasterList:
     rasterObj = arcpy.Raster(raster)
     print raster

     bands = arcpy.GetRasterProperties_management(Mask ,raster, "MEAN")
     print "MEAN VALUE: %s" %bands

Error Msg:

Property type: The input is not within the defined domain.
The value is not a member of MAXIMUM | MINIMUM | MEAN | STD | TOP | LEFT | RIGHT | BOTTOM | CELLSIZEX | CELLSIZEY | VALUETYPE | COLUMNCOUNT | ROWCOUNT | BANDCOUNT | ALLNODATA | ANYNODATA | SENSORNAME | PRODUCTNAME | ACQUSITIONDATE | SOURCETYPE | CLOUDCOVER | SUNAZIMUTH | SUNELEVATION | SENSORAZIMUTH | SENSORELEVATION | OFFNADIR | WAVELENGTH.
Failed to execute (GetRasterProperties).


texts: Remove layer file from MXDs with multiple data frames using ArcPy?

texts:
I have a bunch of mxd's and I want to remove layer file from them. Each mxd has 2 and more data frame. I don't know why it does not remove.

I'm trying this code:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project"
for mxd in arcpy.ListFiles("*.mxd"):
    print mxd
    mapdoc = arcpy.mapping.MapDocument(r"C:\Project\\" + mxd)
    df = arcpy.mapping.ListDataFrames(mapdoc, "Layers")[0]
    df2 = arcpy.mapping.ListDataFrames(mapdoc, "Layers")[1]
    Layer1 = arcpy.mapping.Layer(r"F:\GIS\practical_shapefile\residence\residence.lyr")
    arcpy.mapping.RemoveLayer(df, Layer1)
    print 'remove'
    arcpy.mapping.RemoveLayer(df2, Layer1)
    print 'remove'
    mapdoc.save()
del mxd

This prints out:

Project.mxd
remove
remove
Project2.mxd
remove
remove
Project3.mxd
remove
remove


texts: Is there any reason NOT to use MODIS sinusoidal projection?

texts:
I have a large number of MODIS snow rasters which I am cropping with the extent of a shapefile. Currently, I am using spTransform so that my shapefile has the same CRS as the rasters.

This code works perfectly, it returns a new raster with the extent of the shapefile. My question is: is there any reason to avoid this sinusoidal projection? Mathematically, it should preserve area better than another projection. Are there any disadvantages?

Here is my code, for reference.

prj.aea="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"

##Import shapefile with WGS84 projection
dsn=getwd() # assumes current dir is folder with shapefile(s) to import
ut.mgmt1=readOGR(dsn=dsn,layer=paste(GEO,sep=""))
GEO=spTransform(ut.mgmt1,CRS=CRS(prj.aea))

temp = list.files(pattern="*.tif",recursive=FALSE)
####THIS WRITES A NEW, SMALLER RASTER####
for(x in 1:length(temp)){
    r1 = raster(temp[x])
    test = mask(r1,GEO)
    writeRaster(test,file=temp[x],overwrite=TRUE)
}


texts: ArcGIS JS API Measurement widget change measurement units text?

texts:
The ArcGIS Server JavaScript API has a Measurement widget, which reports the area or length of a shape drawn on the map.

The default metric labels use the American spelling of "Sq Kilometers" and "Kilometers", etc.

var measurement = new Measurement({
        map: map,
        defaultAreaUnit: units.SQUARE_KILOMETERS,
        defaultLengthUnit: units.KILOMETERS
}, dom.byId("measurementDiv"));

measurement.startup();

Is it possible to change the default text "Square Kilometers" shown by default to a custom text (like "km2")?

The Measurement widget doesn't seem to have any customisation options which would accomplish this.


texts: How to Iterate ModelBuilder for Zonal Statistics?

texts:
I want to use ModelBuilder to split a DEM into multiple tiles, then run the iterator on each tile made for the Zonal Statistics tool to find the highest point(s) in each tile and give these as an output.

Each time I try this ModelBuilder crashes.

Any suggestions/ ideas on how to make my model run?


texts: ogr2ogr PostGIS to GeoJSON using batch file

texts:
I am trying to export a PostGIS table to GeoJSON using ogr2ogr.
The following works perfectly when I launch it in Windows command :
ogr2ogr -f "GeoJSON" F:\output.json PG:"host=**** port=**** dbname=**** user=**** password=****" --optfile F:\query.sql`

However, when I save this command as a .bat file and then try to launch it in Windows command, I receive an error :
FAILURE: Unable to open datasource `PG:host=**** port=**** dbname=**** user=**** password=****' with the following drivers.

Is there any way I can make this command work as a .bat file?


texts: Can I open ArcView 3 project files (*.APR) with Qgis? Is there a converter?

texts: I would like to open files from an old ArcView 3.x project with Quantum GIS to work with the data. But this project data is saved as an APR-file, which I can not open with Quantum GIS. Is there any possibility to convert the data (for example to QGS files) to work with it?


texts: Select features matching the values in an array (QGIS 3.10 expression)

texts:
I want to use an expression to select features matching the values of an array.

Select features from a list is done by the IN operator:

"field" IN ('value1','value2')

But how to select features from an array like ['value1','value2']? Do I have to convert it to a list (somehow) or is there a "direct" method of doing so?

QGIS 3.10


texts: Creating DEM after converting MapInfo tab files to ArcGIS Desktop?

texts: How do I create DEM for contours after converting MapInfo tab files to ArcGIS and how do I convert toposheet tab file to ArcGIS?


texts:
I would like to calculate the NDVI for an area to identify vegetation cover over a 32 year period using. However, I am having challenges with the calculations I have attempted with Raster Calculator. I am using images acquired from the USGS Earth Explorer for the years 1988, 2000, 2017, 2018 and 2020. I have created the Multi Band Composites for each year and used the Bands from those Composite images to perform the calculation for NDVI in Raster Calculator.
For the first image 1988 (Landsat 4) I used Band 4 (NIR) and Band 3 (Red), using the equation
Float("1988\compo_1988 - Band_4" - "1988\compo_1988 - Band_3") / Float("1988\compo_1988 - Band_4" + "1988\compo_1988 - Band_3")
However, the results (seen in the screenshot below) were not what I was expecting, I.E values ranging from -1 - 1.

Is there a way I can get the values that are usually associated with NDVI? If not, how can I use these results to identify which areas are vegetation or not?


texts: Split lines by a polygon - Why are "internal" lines split in parts using sf/st_split

texts:
The context: I have a border and routes followed by vehicles. I want to

decompose the routes, cutting them everytime they cross a border, and
for each part, know what are the distances performed.

Simple reproducible example, part 1: input data: as input, one square as polygon for the border, three routes as linestrings:
library(sf)
library(lwgeom)

# Defining the "border", a polygon
four_corners <- matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0),
                       ncol=2, byrow=TRUE)
pts <- list(four_corners)
my_polygon <- st_polygon(pts)
plot(my_polygon)

# Defining the first "route", a linestring
three_points <- matrix(c(-1, 0, 0.5, 0.5, 2, 0),
                       ncol=2, byrow=TRUE)
my_linestring1 <- st_linestring(three_points)
print(my_linestring1)
plot(my_linestring1, add=TRUE)

# Defining the second "route", a linestring
two_points <- matrix(c(0.5, -0.5, 0.5, 0.25),
                       ncol=2, byrow=TRUE)
my_linestring2 <- st_linestring(two_points)
print(my_linestring2)
plot(my_linestring2, add=TRUE)

# Defining the third "route", a linestring
four_points <- matrix(c(0.25, 0.75, 0.25, 0.6, 0.75, 0.6, 0.75, 0.75),
                     ncol=2, byrow=TRUE)
my_linestring3 <- st_linestring(four_points)
print(my_linestring3)
plot(my_linestring3, add=TRUE)

# Group the three routes in a spatial dataframe with ID
g <- st_sfc(my_linestring1, my_linestring2, my_linestring3)
df_routes <- data.frame(id_column = 1:3)
st_geometry(df_routes) <- g

df_routes looks like this:
  id_column                       geometry
1         1 LINESTRING (-1 0, 0.5 0.5, ...
2         2 LINESTRING (0.5 -0.5, 0.5 0...
3         3 LINESTRING (0.25 0.75, 0.25...

With one internal route, one route crossing the border once, and one route crossing the border twice:

Simple reproducible example, part 2: the solution
# Split the route in different subroutes cut by the border
splitted_routes <- st_split(df_routes, my_polygon)

# Decompose the GEOMETRYCOLLECTION of LINESTRINGs in a list of LINESTRINGs
splitted_routes_as_single_linestring <- st_collection_extract(splitted_routes, "LINESTRING")

# Get the length of the LINESTRING
list_of_length = st_length(splitted_routes_as_single_linestring)

# Remove the geometry, in order to get a dataframe as output
splitted_routes_as_single_linestring$geometry <- NULL

# Just merge the length in a column distance in the original dataframe
ID_and_length <- cbind(splitted_routes_as_single_linestring, distance = list_of_length)

This approach works on this simple example and generates the distances of each part on each side of the border:
    id_column distance
1           1 1.054093
1.1         1 1.054093
1.2         1 1.054093
2           2 0.500000
2.1         2 0.250000
3           3 0.800000

In particular, the 3rd route is internal and the output is simply one distance corresponding to the complete distance of the original route.
Now, when using this approach on a real life example, some internal routes are decomposed in small pieces. They should be decomposed as one, but somehow, my approach is not working. When I am extracting LINESTRINGS out of the GEOMETRYCOLLECTIONS coming from st_split, I got several objects, which was not my goal/expectation.
What's happening?


texts: Writing SQL Query to create spatial view?

texts:
I use this sql query to draw a rectangle from a point layer (overview):

create or replace view view_test as
select name, width, height, ST_SetSRID(ST_GeomFromText('POLYGON (('
 ||ST_X(geom)-(width * (scale * 0.1))/2||' '||ST_Y(geom)-(height*(scale*0.1))/2||','
 ||ST_X(geom)-(width * (scale * 0.1))/2||' '||ST_Y(geom)+(height*(scale*0.1))/2||','
 ||ST_X(geom)+(width * (scale * 0.1))/2||' '||ST_Y(geom)+(height*(scale*0.1))/2||','
 ||ST_X(geom)+(width * (scale * 0.1))/2||' '||ST_Y(geom)-(height*(scale*0.1))/2||','
 ||ST_X(geom)-(width * (scale * 0.1))/2||' '||ST_Y(geom)-(height*(scale*0.1))/2||'))'),32632) as geom
from overview;

But QGIS isn't loading the results:

2017-07-14T21:38:12 1   dbname='XXXX' host=XX.XX.XXX.XXX port=5432 user='postgres' sslmode=disable key='width' table="public"."view_test" (geom) sql= is an invalid layer - not loaded

Whats wrong with my code?


texts: Changing ArcMap Display Order of Elements in Same Layer File?

texts: I have a bunch of polygons in a layer file and I want to reorder those polygons so one of them is "sent to back". I am using ArcMap 10.5.


texts: Showing heading on a rectangular map projection

texts:
I'm creating an app which displays a map in the form of a basic lat/long grid along with the user's location and heading. I'm currently plotting the grid lines and user location with an equirectangular projection (i.e. simply treating longitude/latitude as x/y coordinates on a rectangular grid) which works fine, but I have a few questions when it comes to plotting the user's heading.
To calculate the current heading, I'm using the formula under the 'Bearing' section of this page. Whenever I get a location update containing a new latitude/longitude, I put it into the formula along with the previous latitude/longitude to find the bearing between them.
However, when I plot my calculated heading arrow on the map, the arrow does not align with the line drawn between the user's current location and their last location, unless you're at the equator. I'm guessing this is because the equirectangular projection distorts a circle into an ellipse as seen here (i.e. it's not conformal).
So, I'm wondering how I might 'transform' the heading value calculated using the above method so that it accounts for this distortion, and the arrow points in the correct direction? I know that one solution would be simply to do atan2(deltaLat, deltaLon) but I'd much prefer to start with a proper heading value and then transform it to a displayable angle depending on which map projection I choose.
Note that the user is not travelling large distances between each location update, so I wouldn't expect to have to worry about the heading changing along a great circle.
Side note: am I right in thinking that if I used a Mercator projection, I could just plot the heading as-is because the distortion does not occur at extreme latitudes?


texts:
I'm new on PostGIS
I have this

The orange ones are buffers, the red dots are the centroids of each brown polygon
So I've tried to cover the brown polygons by making buffers of the centroids (red dots) and the results are this

How can I get the pink buffers that cover the most parts of the red dots?.
I made something on QGIS picking it by myself manually.

I'm looking for something like this, but no idea how to start the sentence.


texts: Soilgrids layers should but don't all have the same CRS

texts:
I am working with SoilGrids layers downloaded from their file repository (https://files.isric.org/soilgrids/latest/data/). I was trying to stack the layers, but found that, in addition to different extents (see post), layers also use two different coordinate reference systems (CRS). I am using R to process the raster files.
While most layers give:
> library(raster)
> crs(raster("sand_0-5cm_mean.tif"))
CRS arguments:
 +proj=igh +ellps=WGS84 +units=m +no_defs

Toward the beginning of the .vrt file I used to build the geoTIFF, I see
<SRS>PROJCS["Interrupted_Goode_Homolosine",GEOGCS["GCS_unnamed ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Interrupted_Goode_Homolosine"],UNIT["Meter",1]]</SRS>

That CRS is used for most layers I used. However, a few (so far nitrogen_0-5cm_mean.tif, phh2o_0-5cm_mean.tif and soc_0-5cm_mean.tif) use another crs:
> crs(raster("nitrogen_0-5cm_mean.tif"))
CRS arguments:
 +proj=goode +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs

Again, the .vrt files have
<SRS>PROJCS["World_Goode_Homolosine_Land",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["Degree",0.0174532925199433,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","4326"]],PROJECTION["Goode_Homolosine",AUTHORITY["Esri","43059"]],PARAMETER["False_Easting",0.0,AUTHORITY["Esri","100001"]],PARAMETER["False_Northing",0.0,AUTHORITY["Esri","100002"]],PARAMETER["Central_Meridian",0.0,AUTHORITY["Esri","100010"]],PARAMETER["Option",1.0,AUTHORITY["Esri","100036"]],UNIT["Meter",1.0,AUTHORITY["EPSG","9001"]],AUTHORITY["Esri","54052"]]</SRS>

When I plot the maps however, the projections seem identical. From what I found in this PDF from PROJ (proj.org seems down at the moment), +proj=goode is used for the uninterrupted Goode Homolosine projection. However, SoilGrids layers are supposed to all use Interrupted Goode Homolosine which would be +proj=igh. For now, I am just changing the CRS with:
> nitrogen0_5 <- raster("nitrogen_0-5cm_mean.tif")
> sand0_5 <- raster("sand_0-5cm_mean.tif")
> crs(nitrogen0_5) <- crs(sand0_5)
> writeRaster(nitrogen0_5, "nitrogen_0-5cm_mean_fixedcrs.tif")

That seems to work, but I would rather only change the crs if I am really sure about which one it is. Is it on purpose that the crs differ between raster files or is that an error?


texts: Suppress QgsTask "Task Complete" notification

texts:
I'm following the Extending QgsTask example to create a custom QgsTask.

However, I don't want it to display the Task Complete notification* every time the task completes successfully.

Is it possible (preferably using PyQGIS) to suppress the "Task Complete"  notification?

* In Windows 10 at least, the notification appears in the bottom right of the screen above the time (with a ding!).


texts: Add Path to JPGs in Attribute Table based on Unique ID

texts:
I've been tasked with updating attribute table fields with file paths to image files in a folder based on a unique ID. I'm having trouble thinking up a script that will do this? I'm working on a standalone script in ArcGIS Pro 2.3.1 Arcpy.

This is an example of how the table should look. Every JPG has a unique ID that matches the Post_Number field unique ID. Currently, the paths to the JPGs are entered manually into the table, and there are 1000s of them.

I don't think a Calculate Field will work because not all JPGs are stored in the same directory. All I've done so far is create a cursor that searches the Post_Number and Photo fields.

#for each row print the fields
    with arcpy.da.SearchCursor(fc, fields) as cursor:
        for row in cursor:
            print(u"{0}, {1}".format(row[0], row[1]))


texts: Weighted Average of Band Based on another Band in Image Collection [Google Earth Engine]

texts:
Using monthly VIIRS data, I'm trying to estimate a weighted average of radiance where I want to  weight by the number of cloud free days. Each image in the image collection has two bands: "avg_rad" and "cf_cvg." The below code calculates the average radiance in 2015; however, I'm wondering if there's a convenient way to calculate a weighted average, where "cf_cvg" is used as the weight.

var viirs = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG');
var viirs_2015 = ee.ImageCollection(
    viirs.filterDate("2015-01-01","2015-12-31")
);

var viirs_2015_mean = viirs_2015.reduce(ee.Reducer.mean());
var viirs_2015_mean_rad = viirs_2015_mean.select("avg_rad_mean")


texts: Problem selecting features with % operator

texts:
I am having problem to select contours with the following expression:
"ELEV" % 0.05=0
As I understand this should select any contour that the remainder of the division by .05 is 0.
For instance, contours with an elev value of 2.1 or 2.15 or 2.2 or 2.25 should be selected. Nevertheless they are not selected. This affects the rule based symbology or labeling and select features by attributes.
Contours were generated with the QGIS contour extraction tool.
What I am doing wrong?


texts: Problem processing OSM file with osm2pgsql

texts:
I have an OSM file and when I want to import it in my PostGIS base using osm2pgsql, I get the following error:

Processing: Node(130k 8.7k/s) Way(0k 0.00k/s) Relation(0 0.00/s)Assertion failed: (xlon), function StartElement, file parse-xml2.c, line 101.

How can I fix this problem?

I use the following commands that works with other osm files excepting this one:

osm2pgsql "/Users/me/Desktop/my_osm_data.osm" -d my_osm_data -S "/Users/me/Documents/SIG/QGIS/osm2pgsql Styles/TC.style" -W -U postgres --slim -l

Here is the osm file: https://www.box.com/s/5a2qck2fc3tsnszimc4d


texts: Extract areas of multi-part polygons (SpatialPolygonsDataframe) R

texts: I have a SpatialPolygonsDataframe with multiple polygons, some of these have multiple parts (blue shapes below). In GIS-land it would be easy to get the areas of the individual parts of a multi-part polygon but it seems difficult to do this in R?


texts:
I have this code (below) for plotting some information in a map, but it doesn't work, in the map appear only the error "sppanel", why I solve that?

library(GWmodel)

library(RColorBrewer)

data(DubVoter)

gwss.1 <- gwss(Dub.voter,vars = c("GenEl2004", "LARent", "Unempl"),
               kernel="bisquare", adaptive=TRUE, bw=48)

gwss.mc <- montecarlo.gwss(Dub.voter,vars = c("GenEl2004", "LARent",
                                              "Unempl"), kernel="bisquare", adaptive=TRUE, bw=48)

gwss.mc.data <- data.frame(gwss.mc)

gwss.mc.out.1 <-ifelse(gwss.mc.data$Corr_GenEl2004.LARent < 0.975 &
                         gwss.mc.data$Corr_GenEl2004.LARent > 0.025 , 0, 1)

gwss.mc.out.2 <-ifelse(gwss.mc.data$Corr_LARent.Unempl < 0.975 &
                         gwss.mc.data$Corr_LARent.Unempl > 0.025 , 0, 1)

gwss.mc.out <- data.frame(Dub.voter$X, Dub.voter$Y, gwss.mc.out.1,
                          gwss.mc.out.2)

gwss.mc.out.1.sig <- subset(gwss.mc.out, gwss.mc.out.1==1, select =
                              c(Dub.voter.X, Dub.voter.Y, gwss.mc.out.1))

gwss.mc.out.2.sig <- subset(gwss.mc.out, gwss.mc.out.2==1, select =
                              c(Dub.voter.X, Dub.voter.Y, gwss.mc.out.2))

pts.1 <- list("sp.points", cbind(gwss.mc.out.1.sig[,1],
                                 gwss.mc.out.1.sig[,2]), cex=2, pch="+", col="black")

pts.2 <- list("sp.points", cbind(gwss.mc.out.2.sig[,1],
                                 gwss.mc.out.2.sig[,2]), cex=2, pch="+", col="black")

mypalette.gwss.1 <-brewer.pal(5,"Blues")
mypalette.gwss.2 <-brewer.pal(6,"Greens")

map.na <- list("SpatialPolygonsRescale", layout.north.arrow(), offset = c(329000,261500), scale = 4000, col=1)

map.scale.1 <- list("SpatialPolygonsRescale", layout.scale.bar(),
                    offset = c(326500,217000), scale = 5000, col=1, fill =
                      c("transparent", "green"))

map.scale.2 <- list("sp.text", c(326500,217900), "0", cex=0.9, col=1)

map.scale.3 <- list("sp.text", c(331500,217900),"5km", cex=0.9,col=1)

map.layout.1 <-list(map.na,map.scale.1,map.scale.2,map.scale.3,pts.1)
map.layout.2 <- list(map.na,map.scale.1,map.scale.2,map.scale.3,pts.2)

library("sp")
install.packages("lattice")
library("lattice")

X11(width=10,height=12)
spplot(gwss.1$SDF,"Corr_GenEl2004.LARent",key.space = "right",
       col.regions = mypalette.gwss.1,at=c(-1,-0.8,-0.6,-0.4,-0.2,0),
       par.settings = list(fontsize=list(text=15)), main = list(label="GW
correlations: GenEl2004 and LARent", cex=1.25), sub=list(label="+
Results of Monte Carlo test", cex=1.15), sp.layout=map.layout.1)


texts: How to convert stylemap with adduniquevaluerules to openlayers 3?

texts:
I am new to openlayers. I have following code in openlayers 2.

for (var i=0;i<wfslayers.length;i++) {

    var subteStyleMap = new OpenLayers.StyleMap({
        "default": new OpenLayers.Style({
            'strokeWidth': 1,
            'strokeOpacity': 1,
            'fillOpacity': 0.5,
            'fillColor':"orange",
            'strokeColor':"orange",
            'pointRadius': 3
        }),
        "select": new OpenLayers.Style({
            fillColor:"#00FFFF",
            strokeColor:"#00FFFF",
        })
    });

    var object=wfslayers[i];
    var styles=[];

    for (var j=0;j<wfsStyles.length; j++) {
        if(object.wfs_id==wfsStyles[j].wfs_id) {
            styles=wfsStyles[j].data;
            break;
        }
    }
    if (styles[0]) {
        var category=styles[0].attribute_name;
        var lookup={};
        for (var k=0;k<styles.length;k++) {
            var style=styles[k];
            lookup[style.attribute_value]={
                strokeColor: style.strokecolor,
                strokeWidth: style.strokewidth,
                fillColor: style.fillcolor,
                pointRadius: style.pointradius
            };
        }
        subteStyleMap.addUniqueValueRules("default", category, lookup);
    }
    var layer = new OpenLayers.Layer.Vector(object.name, {
        layerid: object.label_text,
        projection: new OpenLayers.Projection("EPSG:"+object.srid),
        styleMap: subteStyleMap,
    }
    wfscatalog.push(layer);
}

With the help of this: https://stackoverflow.com/questions/28093975/openlayers3-make-choropleth-map-with-geojson , I tried to write the equivalent openlayers3 code. But I need help on stylemap. Please help.

for (var i=0;i<wfslayers.length;i++) {

    var styleNormal = new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(255, 128, 64, 1)',
            width: 1
        }),
        fill: new ol.style.Fill({
            color: 'rgba(255, 128, 64, 0.5)'
        })
    });

    var styleSelect = new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: '#00FFFF'
        }),
        fill: new ol.style.Fill({
            color: '#00FFFF'
        })
    });
    var select = new ol.interaction.Select({style:styleSelect});

    var object=wfslayers[i];
    var styles=[];

    for (var j=0;j<wfsStyles.length; j++) {
        if(object.wfs_id==wfsStyles[j].wfs_id) {
            styles=wfsStyles[j].data;
            break;
        }
    }
    if (styles[0]) {
        var category=styles[0].attribute_name;
        var lookup={};
        for (var k=0;k<styles.length;k++) {
            var style=styles[k];
            lookup[style.attribute_value]={[
                new ol.style.Style({
                    stroke: new ol.style.Stroke({
                        color: style.strokecolor,
                        width: style.strokewidth
                    }),
                    fill: new ol.style.Fill({
                        color: style.fillcolor
                    })
                });
            ]};
        }
    }
    var layer = new ol.layer.Vector({
        projection: "EPSG:"+object.srid
        style: function(feature, resolution) {
            return lookup[feature.get(category)];
        }
    });
    layer.layerid=object.label_text;
    wfscatalog.push(layer);
}


texts: How to transfer data and replica information from a file geodatabase to ArcSDE?

texts:
I've got an arcgis 10.1 file geodatabase that contains a number of featureclasses and a number of one-way child replicas.

Is there any tool that I can use to transfer the data and replica information verbatim into an ArcSDE geodatabase?  I'm giving the replica file geodatabase to an external person and I want him to easily be able to move the data and replica information to an ArcSDE database which he will then version and publish as a geodata service that I'll replicate to.

The Extract Data Wizard in the distributed data toolbar does half the job by transferring the featureclasses verbatim (including GlobalIDs).  But replica information does not appear to be transferred.

Create Replica doesn't appear to help in this situation.


texts: Joining a .dbf table to a feature class in a personal geodatabase only matches 90% of the records

texts:
The matching ID fields are identical, there are no extraneous spaces or anything like that. I exported the feature class and table to a gdb, and it joins all but one record (still no difference in the IDs). I need to do this for multiple feature classes. It makes no sense to me that it works in a gdb but not in an mdb.

And yes, the data has to be in a mdb - it's for modelling software.

EDIT:

The field types are identical - string with lengths of 32 in both. No spaces, no dashes, fields start with a letter.

So I just tried joining the exported feature class with the original dbf, and that only matches 90%. But when they are in the same gdb, it does all but one record.

"C_LRC_RCD2_100" is a matching record.

"C_LRC_RCD2_MT_20" is a non-matching record.

Using ArcGIS 10.3.1

Could it have something to do with the feature class join field not being indexed?

It is a simple, one-to-one join

Text from 'Validate Join':

All field and datasource validation tasks were completed successfully.

The number of matching records for the join:
- 288 of 317 records matched by joining [INFO_COND_ID] from  with [ID] from .
Matching records may not appear in table view due to join validation errors.

"C_LRC_RCD2_MT_2" is a matching record, "C_LRC_RCD2_MT_20" is a non-matching record:


texts: Downloading ArcGIS Online data, and editing in ArcGIS as layers, not separate geodatabases?

texts:
My goal: I'm trying to download layers from ArcGIS Online and put them into ArcGIS Desktop (10.3.1) as layers and edit them. One is a feature layer, another is a hosted feature layer. Most of the information I want to use is in the feature layer that is not hosted.

My problems:

The feature layer that isn't hosted cannot be downloaded for use in ArcGIS, is it possible to change it to a hosted feature layer in an editable format (i.e. not a GDB  - see below), OR do I need to start over?

This is what opens when I download the map or the layers:

I have very little experience with ArcGIS Online, though I am familiar with ArcGIS (not python). I am the owner of all the layers and the administrator of the organization.


texts: How to correctly style GEOJSON (ESRI shapefile) for Geoserver?

texts:
I have a geojson file, where features are styled - fill and stroke colours - like this:

{"type":"Feature"
,"geometry":{"type":"Polygon","coordinates":[[[15.50487,50.72638],....,15.50487,50.72638]]]}
,"style":{
      "stroke":"#FF0000",
      "stroke-opacity": 1.0,
      "fill":"#FB7F7F",
      "fill-opacity":0.8
 }
}

I convert the file to shp format using ogr2org like this:

ogr2ogr -f "ESRI Shapefile" as.shp as.geojson -nlt POLYGON -skipfailures

When I create a layer using the shp file, all features are displayed in grey although the style specified in original geojson is shown in table below the Geoserver preview map (see the picture for the detail).

What am I doing wrong?


texts: Plot Canadian (Alberta Province) township, range, and section data on Google Maps

texts: We are working on a project where we need to plot some part of Canada (Alberta) over Google Maps. We can plot maps on Google with coordinates (lat/long) or geometry values. But our client provided us data in form of Township, Meridian, Range, Section, and LSD (legal subdivision or quarter). Could someone out there guide us how can we plot this data on Google Maps? Client doesn't want to use any converter to convert this data to coordinates as he assumes that TRS data has its own curvature, etc.


texts: Python 2.7.13 cannot import ArcPy in ArcGIS Desktop 10.5.1?

texts:
I am using ArcGIS 10.5.1 and Python 2.7.13. I have also ArcGIS Pro on the same computer. Our license is in network server in our company.

Import arcpy is currently not working.

I had some tools to process shapefile and do some calculations and I have also implemented Tkinter GUI. It worked before. Suddenly it stops working, it cannot import ArcPy any more and automatically terminates the program without showing any error message.

RESTART: C:\Users\SPervin\Desktop\Tools_Biotop_FFH_berechnung\Biotop_berechnung_tkinter.py
RESTART: Shell
import arcpy
RESTART: Shell

It is just terminating the the process.

Example screenshot for system paths:


texts: Get only fields with datatype int in PyQGIS

texts:
With the following line of code, I get a list of field names for the layer layer.

 fields = layer.fields().names()

However, I would like to get only those fields with the datatype Integer. How would I do that? I checked the class QgsFields but I don't see a function that relates to the datatype. Is there a way to do that?


texts:
I'm trying to import WKT (http://polygons.openstreetmap.fr/get_wkt.py?id=2413177&params=0 - boundary of one of districts of Łódź) into PostGIS database and convert it to geometry. I've created table 'districts' with 'the_geom geometry', id and name of the district.
When I tried to insert data with

INSERT INTO districts (name, the_geom) VALUES ('Polesie', ST_GeomFromText('<<<WKT goes here>>>', 4326))

the the_geom field is left empty. In fact, even query like SELECT ST_GeomFromText('...', 4326) returns empty query. I also tried to create text field in the table, insert WKT data in there and then run UPDATE districts SET the_geom=ST_GeomFromText('...', 4326) but it doesn't work.
I tried to use query from http://www.postgis.org/docs/ST_GeomFromText.html and it worked.

Could anyone help me with this one? Am I doing something wrong, or the WKT is just too big?

Definition of my table:

CREATE TABLE dzielnice
(
  name character varying(30),
  the_geom geometry,
  wkt text,
  id bigserial NOT NULL,
  CONSTRAINT dzielnice_pkey PRIMARY KEY (id),
  CONSTRAINT dzielnice_id_key UNIQUE (id)
)


texts: Display name when hover over plugin icon in QGIS

texts:
I have made a QGIS plugin (Python) and when I hover my mouse/cursor over the icon in QGIS the name of the plugin seem to be "Plugins Toolbar". I have searched all .py scripts in the directory of the plugin and can't find any reference to "Plugins Toolbar".

How do I change the name of this "icon text"?


texts: Opening feature class created in_memory using ArcObjects?

texts:
How do you open an IN_MEMORY featureclass? The code below is not working:

            IGPUtilities2 gpUtils = new GPUtilitiesClass();
            IFeatureClass featureClass = gpUtils.OpenFeatureClassFromString(@"in_memory\MergeHT");
            IFields2 fields = featureClass.Fields as IFields2;
            IField2 nameField = null;
            for (int fieldIndex = fields.FieldCount - 1; fieldIndex > -1; fieldIndex = fieldIndex - 1)
            {
                nameField = fields.Field[fieldIndex] as IField2;
                if (nameField.Name != "FID" && nameField.Name != "Shape" && nameField.Name != "OBJECTID" && nameField.Name != "SUBTYPECD")
                {
                    featureClass.DeleteField(nameField);
                }
            }

It gives below error:


texts:
I have created a vector layer by reading a geojson file and I have a SLD file created for this layer. I am wondering how do I link the SLD file to the layer. I had a read of this example on Openlayers website:

http://openlayers.org/dev/examples/sld.html

I am having a bit of trouble understanding the code for this example- which are using radiobuttons to change the style to the selected button. I believe the linking in my situation will be much simple than what is showing in the example- just linking one SLD file to one single vector layer.

Any guidance will be greatly appreciated.


texts: GeoExplorer dependencies are missing after git clone

texts:
When I set up GeoExplorer in development mode today I had to manually clone the external js libs that are referenced in git (buildkit, stick, geoext, gxp, openlayers, printpreview, proj4js).

They were not included after running "git clone git://github.com/opengeo/suite.git", which prevented me from building GeoExplorer with maven. After adding them manually everything works. What am I missing in my git/maven configuration to skip this step?


texts: Exporting to DXF using Python doesn't give layername compared to using QGIS File export

texts:
I'm trying to export a DXF file with python in QGIS, but I get different results when I use "Project->Import/export->Export project to DXF...-> and then export to a DXF file"
Compared to my Python script, the export works fine and I get a DXF file which contains all information except the layer name.
When I place the DXF file in my CAD program the one I make with python does not contain the layer name. While the one I export using "Project->Import/export->Export project to DXF...-> and then export to a DXF file" does.

Below is the python code i'm using
Does anyone know if this is possible using Python?
output_name ="ABCD"
output_filename = os.path.join(directory + dxfmap, output_name + '.' + "dxf")

options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "DXF"
options.fileEncoding = "CP1252"
options.onlySelectedFeatures = False
options.destCRS = layer2.crs()
options.layerName = "ABCD Layername"
options.skipAttributeCreation=True
options.setLayerTitleAsName=True

writer = QgsVectorFileWriter.writeAsVectorFormat(
        layer2,
        output_filename,
        options
)


texts: Running Lookup tool in ModelBuilder using field that was joined in previous step

texts:
I'm trying to generate a model (ModelBuilder in ArcMap) that calculates the Landslide Susceptibility Index with standardized inputs (Conditioning factors in raster format and dbf's with weights). The last step is the raster calculator where I basically use a weighted sum expression
(Slope * slope_weight + Geology * geology_weight + Aspect * aspect_weight + ...)

For this I need the several rasters with the class weights (each class has a different weight) and I used the Join field  tool (Input and Output join field=Value and field to be joined=weights) and the Lookup tool (Lookup field=weights from the previous join) to generate the raster layers with class weights (float) as pixel values.
The issue is that the join field is not recognized by the Lookup tool and therefore the model doesn't run. I supposed that if I set the Join field process as a precondition of the Lookup tool, the model would start running and the lookup tool would "wait" for its precondition (joined field), but that is not the case :(.
Is there anything I can do to set (in ModelBuilder) make this model run as I planned it?


texts:
I'd like to identify the layers that are currently active (which basemap is currently in view and which overlays are currently turned on).

Was hoping to find something akin to getActiveLayers in the Leaflet API (http://leafletjs.com/reference.html/) but I'm not seeing anything...


texts: Maturity of Mapserver mapscript C# bindings?

texts:
I've been thinking about possibly porting an internal custom GIS application built on top of Python / PyQt / Mapscript / GDAL-OGR / PostGis to C#.

I was wondering about the maturity of the C# mapscript bindings (and to a lesser extent, the GDAL/OGR bindings). Can they be considered as mature as the current (very nice) Python bindings since I assume both are SWIG generated?

(The most time intensive part will probably be redoing the current PyQt GUI parts)

Thanks in advance for any insight / experiences concerning this topic!


texts: How to eliminate small features in a grid?

texts:
(Update) My final goal is to take some logged point data (1 point about every 20 ft or so there are over 15,000 points) and interpolate that point data.  Then I would like to convert the interpolated raster into distinct polygon zones based on the number of classes the data is in using the Natural Breaks classifications.  So if I want it in 5 classes, there will be several polygons that make up these 5 classes, but within the tabular data I can see which of the 5 classes each polygon belongs to.  In addition to this, I would like to eliminate the smaller polygons created by the interpolation, like say anything under 1 acre in size.  Maybe there is a GP process I can run to do this?  Let me know if I need to explain this better, and thanks for the help so far.

I was thinking the "output cell size" setting would be the trick here but it defaults to a very small number (0.00003) and if I make it any bigger like even say "0.003" it results in the interpolated raster having virtually no data where both high and low values is a really high number and very low number.

If I use the default settings I get a decent kriged raster but I would like to adjust to exclude smaller areas, like say less than 1 acre in size.  Basically I want the kriging interpolation to ignore areas that are smaller than an acre or even a 1/2 acre.  Currently it is showing areas that are easily less then 0.001 acres in size of which I would like it to ignore those. Also if I zoom into it looks like the cell size grids are 8ft by 8ft.  How do I get that larger to say 100 ft by 100 ft?  Like I said I have tired putting larger number in the output cell size but get back basically nothing but extreme values.  I must be doing something wrong.  Any help would be appreciated.


texts: This is something I notice recently. I'm at a stage where I'm finally moving on from ArcGIS 9.3.1 to 10.1. I noticed an assiue during this transition. In order to view my metadata, I will need to upgrade my FGDC metadata to the new format which i am aware of (ArcGIS Metadata format). However, once I upgraded, I noticed a clear issue in which some of the content (especially time and date got modified).You can see this in the xml file. Any idea how i could work around this issue ? Ideally i'll like upgrade so that users can read the metadata in ArcCatalog and at the sametime preserve all the content within my FGDC metadata. If not possible, what's the work around solution to this ?


texts: Conversion From Transverse Mercator To Longitude & Lattitude

texts:
I have to transform (7-parameters datum transformation) coordinates from northings & eastings to lattitude & longitude.

Following are details of source coordinate system(i.e. northings & eastings one) :

1) Projection Details:
 1.1)Projection Name
 1.2)Origin Lat/Lon
 1.3)Origin North/East
 1.4)Scaling Factor
2)Datum Name: AGD84
3)7-parameter Datum Details:
 3.1)Rotn X:0.292
 3.2)Rotn Y:-0.443
 3.3)Rotn Z:-0.277
 3.4) Trans X:117.763
 3.5) Trans Y:-51.51
 3.6)Trans Z:139.061
 3.7)Scaling Factor:-0.03939657799319541

Please note that I am using "Proj4Js" for the above mentioned transformation for which following is the "defs" settings:

***Proj4js.defs["<EPSG-CODE>"] = "+proj=<someProjectionImplementation>
+towgs84=-117.763,-51.51,139.061,0.292,-0.443,-0.277,-0.03939657799319541 +no_defs";***

Could anyone confirm that is there any issue in this setting?

Thanks,
Gajendra


texts: Using OpenLayers for projection of static image on worldmap?

texts:
Using the latest OpenLayers, I'm trying to overlay a simple static image over the regular OSM streetmap. Although there are examples how to overlay an simple static image which is in the same projection as the map, or images that use a predefined projection, I have a hard time finding an example with a custom image.

So, what do I have; a normal jpg image, with a certain width & height, and 3 reference points. For example, a reference point will give me that pixel [10, 30](top,left) is at geocoordinate [-73.9984252,40.7197263](EPSG:4326). In theory using the 3 points it should be possible to get a projection that transforms the image so it overlays the map.

Looking at the supported options of openlayers, where for example rotation of the static image already does not seem be supported, I think I need to generate my own proj4 projection, but I have no clue how to start with this. I did manage to reposition the image using lat_0 & lon_0, but that's about it.

Am I on the right track?

Are there better ways to approach this?

Can you give me a hand with the projection?

I've created a codepen.io with an example. Note that my usecase is on level building floorplans, but the idea is the same.


texts: Can I add Google Maps as WMS layer to GeoMedia?

texts:
My primary target is to make a WMS layer for Google Maps and load it to GeoMedia professional.

I read a lot for the past few months and I can't do it only with reading, so I decided to install GeoServer and give it a try. It didn't work so far.

I'm loading a shapefile which I generated from GeoMedia, but I can't figure how to add a map (Google Maps) to my WMS.

It will be nice if someone can help me with some information or example how to do it (if it is possible). I'm 90% sure that it is possible, but I don't know how to do it.

First I created a shapefile from GeoMedia. After that I loaded the same shapefile to GeoServer, I added some styles and I created my WMS.  After that I loaded the WMS to GeoMedia.

The problem is that I want to have a Google Maps or Bing maps as a background for my WMS.

I created an HTML file of Google Maps and added my WMS on it, I put it in my geoserver/www dir and started it from my browser as "localhost:8080/geoserver/www/myfilename.html" and its working fine.

The problem is when I tried to add the same file to GeoMedia - there is an error. I clicked on more to see the problem and I get this:

  Version negotiation failed. The WMS does not support the versions
  supported by the data server (1.3.0/1.1.1/1.1.0/1.0.0).

I made a screenshot of the error message:

.

Can someone tell me why I receive this error?


texts: Display geojson on WMS images in openlayers 3


texts: Align content in an attribute table to the bottom of the frame

texts:
I have an attribute table in a page in one frame with a certain dimension (width/height). The content of this table is not enough to fill the height of the frame and is aligned to the top. Is any of this possible?

1) keeping the fixed dimensions of the frame, align the content of the table to the bottom of the frame?

2) having the frame height automatically calculated based on the content? (i.e. based on the number of rows?)

The table is filtered using the expression "ref_pln" = attribute(@atlas_feature,'ref_pln'), so the table displays rows if the value of their "ref_pln" attribute is equal to the value of the "ref_pln" being evaluated by the atlas.

How can I fill "table count" in the atlas using this filter?


texts: Read the band value and Write the band value in Raster(.tif) using QGIS python

texts:
I have written a plugin in QGIS 3 to read the band value of a pixel or from a textbox and update the value to the selected pixels. When I'm updating a single pixel or smaller group, the script is working fine. When I'm updating for a large group of pixels, the script terminating from the half way and shows the error like  ".tif file is not regonised file format..."

In the below screen, I have entered band value 5 and updated the same for small area. When I tried the same for bigger area, it applies partially and shows the error as mentioned above.

Below is the code used for read, write and close the raster (.tif) file.

def change_cell_value(self, vals, x=None, y=None):

    if not x:
        x = self.px
    else:
        pass
    if not y:
        y =self.py
    else:
        pass
    self.gdal_raster = gdal.Open(self.raster.source(), GA_Update) #Read
    for nr in range(1, min(4, self.bandCount + 1)):
        s = vals[nr-1]
        rbuf = s.tostring()
        self.bands[nr]['data'].WriteRaster(x, y, 1, 1, rbuf, buf_type=self.bands[nr]['qtype']) #Write
        self.bands[nr]['data'] = None
        self.bands[nr] = None
    self.gdal_raster = None  #Close
    self.raster.setCacheImage(None)
    self.raster.triggerRepaint()
    self.prepare_gdal_raster(True)


texts: GeoTools CRS.findTransform suddenly failing

texts:
I have a body of users using an Eclipse RCP-based Java desktop app  (32-bit JRE1.7, Windows 7) that uses GeoTools for world/screen projection handling.

This morning the app has started failing for just one user. It worked fine on Friday. Significantly, CRS is not finding the transform, so _degs2Metres takes a null value:

_worldCoords = CRS.decode("EPSG:3395");
final CoordinateReferenceSystem worldDegs = CRS.decode("EPSG:4326");
_degs2metres = CRS.findMathTransform(worldDegs, _worldCoords);

This happens for just one user, on a network of users with (apparently) identical PCs.  He tried restarting his PC, but it didn't fix it.

He has prior versions of the app running back to Oct 2018, and all of them now fail in the same way.  We also tried opening 10-year old datafiles, and they won't open.  This makes it look like the problem is outside this app.

It really looks like something we're unaware of has changed with his PC config.  He's not aware of anything changing over the weekend, nor are his IT support staff.

Has anybody encountered a class-path, or some MS Windows config change that has caused the CRS.findMathTransform() method to fail?

I believe GeoTools uses an H2 database to store the EPSG codes, so I guess his PC could have some other H2 database driver that is being used instead of the GeoTools one.

Any other suggestions?


texts: Can't get the where clause in Select layer by attribute working!

texts:
I am trying to develop a script that searches for a particular ID in a user selected field. To this end I have used the select layer by attributes tool to search for a derived number in a field specified by the user.I cannot, however, after numerous attempts, get the where clause right. I have posted my latest attempt below:

import arcpy
from arcpy import env

work_space = "L:\\DRAINAGE\\D_GROUP6\\STAFF\\CJE\\EDRNb\\EDRN_holding.gdb"

env.workspace = work_space

arcpy.env.overwriteOutput = True

selected_li = "L:\\DRAINAGE\\D_GROUP6\\STAFF\\CJE\\EDRN_holding.gdb\\EDRN7\\EDRN_LINK"
OB_ID = "OBJECTID_1"

rows = arcpy.SearchCursor(selected_li)
for row in rows:
    id_LCC = row.getValue(OB_ID)
    id_TNDE = row.getValue("TO_DRNNODE")
    if id_TNDE == " " or id_TNDE == None:
        arcpy.MakeFeatureLayer_management(selected_li, "feature_layer10")
        arcpy.SelectLayerByAttribute_management ("feature_layer10", "NEW_SELECTION", '"' + str(OB_ID) + '" =' + "'" + str(id_LCC) + "'")

So, if anyone has any ideas please do let me know.


texts: Is it possible to call FME transformers (or workbenches) from my python script?

texts:
I have a script mostly using arcpy, and in the middle of this script I convert a load of polylines to polygons:

arcpy.ET_GPPolylineToPolygon("Contours_L_1m.shp", "Contours_G_1m.shp", "FALSE")

This worked when it was a very small data set, but with lots of data it takes ages then crashes because it runs out of ram.

Fortunately I found that using the AreaBuilder transformer in FME will do the same thing really quickly and without crashing. The problem here is that I now have to run the first half of my script, stop, convert my contours in FME, then run the second half of my script. So what I would like would be (ideally) to be able to run the AreaBuilder transformer from the ArcGIS python window, or (failing that) run the workbench from the python window, wait for it to finish, then continue the script.

Is this at all possible? and how would I go about achieving it?


texts: How to combine multiple overlapping polygon shapefiles into 1 - retain top most layer

texts:
Working with ArcGIS 10.2...

I'm looking at geological data and have 4 separate polygon shapefiles. They all fall ontop of one another like a stack.

I need to create a single dataset from the four separate shapefiles that only has the visual geology from the display order.  I do not need to include information of additional overlapping features that fall underneath a feature.

I've tried to merge, union, append, update. But I don't think it happening in the correct order.


texts: Using Arcade to symbolize based on edit date AND additional field

texts:
I currently have a script written to symbolize my features based on the date they are edited (below). I want to add a line before this block that evaluates another text field, Mnt_Cdn.  If Mnt_Cdn contains the value "Dead", I want the code to return the value Dead.  If not, I want each record to be evaluated by the date script below.
How would I phrase this?
var days_dif = DateDiff(Now(), water_date, "days");

var result = "water status";
if (days_dif <= 14) {
    result = "Complete";
} else { (days_dif > 14)
    result = "Needs Water" + $feature["plt_year"];
}
return result;


texts:
I cannot offset a line to a specific decimal place, (ie only whole numbers can be inputted)
Am I missing something here??


texts: Displaying intersection from FME LineOnAreaOverlayer transformer in FMEViewer and ArcMap is different?

texts:
I have one very interesting issue which I do not know how to solve/understand. In FME I used the LineOnAreaOverlayer transformer to get new intersected lines, these some are inside and outside polygons. When I visualize it with FMEViewer and visualize them with ArcMAP I get different results.

In ArcMAP it is visible like line was not intersected according to the polygon, but at the same time in FMEViewer everything is fine.
help please!

Both representing same line and area.

And I tried to do intersection within ArcMap and intersections occures at the same place as with FME:


texts:
I have written a script that searches through microstation CADD files, finds ones that conform to a certain naming convention, and pulls them into a geodatabase using the CADtoGeodatabase tool. This script works fine for about an hour and then stops working. No error messages show up in the console. A window just pops up saying that python has stopped working.

I know that there is about 3000 files that need to be converted and my script crashes after about 220. Before trying the conversion I ensure that the file path has no spaces and I have removed all print statements from the script. Adding these two steps allowed the script to run a little longer (originally crashes after about 200 files).

What might be causing this to happen?

import os, logging, re, arcpy

LOG_FILENAME = 'CADcopy.log'
logging.basicConfig(filename=LOG_FILENAME, level=logging.DEBUG)
rootDir = "CADDirectory"

exten = 'namingConventionMatchRegex'
arcpy.env.workspace = "workspace"
counter = 0
def step(ext, dirname, names):
    #print dirname

    x = re.search("xample",dirname.lower())
    if(x):
        return #Do not try to process example data
    for name in names:

        r = re.search(ext, name.lower())
        if(r):
           global counter
           counter = counter + 1
           fileName = os.path.join(dirname, name)
           s = re.search(" ",fileName)
           if(s):
               #print("{} skipped because file path contains space character.".format(name))
               logging.debug("Unable to convert ,{},. File location: ,{}".format(name, os.path.join(dirname, name)))
               return

           #This block gets basic information about the file from its directory tree and incorporates it into the new dataset name
           if name.lower().endswith("endingType1"):
               planType = "Type1"
           else:
               planType = "Type2"

           split = fileName.split('CommonDirectory')
           split2 = split[1].split('\\')
           table = split2[1]+split2[2]+'_'+str(counter)+planType

           #Try to convert to geodatabase. If it already exists, do nothing, else log the filename and path for later review
           try:
               arcpy.CADToGeodatabase_conversion(fileName, arcpy.env.workspace, table, 100)
           except:
               msg = arcpy.GetMessages(2)
               if(msg.startswith("ERROR 000258")):
                  pass
                  #print("File already exists")
               else:
                  #print("Failed to convert {} at filePath {}".format(table, fileName))
                  logging.debug("Failed to convert ,{}, at filePath ,{}".format(table, fileName))

os.path.walk(rootDir, step, exten)