texts: R dismo randomPoints can't handle projected data?

texts:
I'm trying to generate random points from a raster using randomPoints from dismo in R. However, when I overlay the points on top of the raster, it looks like the points aren't in the same projection as the raster:

r <- raster("C:/raster.tiff")
crs(r)

  CRS arguments:
     +proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80
    +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

bg<-randomPoints(r,10000,p=preloc, ext=extent(r),lonlatCorrection = TRUE)

Is there any argument that I am missing to make this line up correctly?


texts: Convert COVID data into data points using ArcGIS Pro

texts:
I am trying to create a script that takes Covid-19 data in a CSV file format and put it into a ArcGIS Pro project using ArcPy. Here is my code. However, I am getting an error saying the parameters are not valid. It says that my in table name does not exist even though it does. What am I doing wrong?
import arcpy

arcpy.env.workspace = r"C:\Users\butle\Downloads\final project"
#gives the location of the csv file
csv_file = r'C:\Users\butle\Downloads\County-Case-Counts.csv'

inTable = "County-Case-Counts.csv"
outLocation = r'C:\Users\butle\Downloads\final project\covid-19data.gdb'

arcpy.AddField_management('Conty-Case-Counts.csv','counties',field_is_nullable="NULLABLE")

expression = arcpy.AddFieldDelimiters("County of Residence") + " = 'covidata' "

arcpy.TableToTable_conversion("County-Case-Counts.csv", r'C:\Users\butle\Downloads\final project\covid-19data.gdb','covid-19data')

inFeatures = "covid-19data"
joinTable = "County of Residence"
joinField = "US_Counties"
expression = "County-Case-Counts.US_Counties = 1"
outFeature = "countyjoin"

county_joined_table = arcpy.Addjoin_management(inFeatures, joinField, joinTable, joinField)

arcpy.SelectLayerByAttribute_management(county_joined_table, "NEW_SECTION",expression)

arcpy.CopyFeatures_management(county_joined_table, outFeature)

arcpy.CalcuateField_management("County-Case_Counts.csv", "countyjoin")


texts:
I am doing versioning of sde. I have done this using the ArcMap. I want to look at the delta tables in sql server.

In sql, I don't see the delta tables. I just see few tables named as a144.

How do I look at the delta tables?


texts:
I am trying to create a mosaic of all the images in file_list, Using:

> vrt_options = gdal.BuildVRTOptions(resampleAlg='cubic', addAlpha=True)
>     gdal.BuildVRT('C:\\my.vrt', file_list, options=vrt_options)

And file_list is a path of all the files i want merged.

 'file_list'
    ['C:\\fasih\\Test_combine\\1974106.mat.tif',
     'C:\\fasih\\Test_combine\\1974107.mat.tif',
     'C:\\fasih\\Test_combine\\1974108.mat.tif',
     'C:\\fasih\\Test_combine\\1974109.mat.tif',
     'C:\\fasih\\Test_combine\\1974110.mat.tif',
     'C:\\fasih\\Test_combine\\1974111.mat.tif',
     'C:\\fasih\\Test_combine\\1974112.mat.tif']

I get this output in Notebook, but nothing is saved in my drive.

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x000002180B46F5D0> >


texts: CORS error when using Google Elevation API

texts:
I'm trying to use Google elevation API with an AJAX request. It doesn't work and when I do console.log(ajax.responseText) I'm getting this error :

Cross-Origin Request Blocking: The Same Origin policy does not allow you to view the remote resource located at https://maps.googleapis.com/maps/api/elevation/json?locations=39.7391536 , -104.9847034 & key = AIzaSyBXSbBWMZtGL2whXVC7T8wFhQcwcKGRuF4. Reason: The CORS header "Access-Control-Allow-Origin" is missing.

Does someone knows how to solve this problem ?


texts: How to find the actual projection name (code) in a PROJCS text data?

texts:
Given a PROJCS data

PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223560493,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",-25.91105555555576],
    PARAMETER["standard_parallel_2",-32.72172222222201],
    PARAMETER["latitude_of_origin",-29.51611111111111],
    PARAMETER["central_meridian",-52.09266666666708],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-725301.460735953995027,572835.914701772970147)
Pixel Size = (169.335263023642995,-169.335263023642995)

How do I find the actual projection name?
I can see EPSG:7030, EPSG:6326, EPSG:4326 and EPSG:9001.
What is the actual EPSG:xxxx of this projection?

Comparing to the known EPSG:900913 I can conclude the last AUTHORITY is naming the Projection:

PROJCS["WGS84 / Google Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Longitude", EAST],
    AXIS["Latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Mercator_1SP"],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH],
  AUTHORITY["EPSG","900913"]]

Is this correct? So my example will be EPSG:9001?

I'm in troble to create a layer in my geoserver using a GeoTIFF with this projection. While I can preview the layer fine, I can't put it to work in a real map or grouping with other layers. It always go to 0,0


texts:
We are having trouble deleting a version made in our 9.3.1 sde.  When we attempt to delete it, ArcCatalog returns a message saying Unable to Delete - Target State not found in States table. Upon looking at the states table there is no reference to that particular version.

Is there any other way to delete this version? We will be upgrading to Version 10.1 soon, will this cause a problem or can it be deleted during the migration process?


texts:
Can you please write me sample loop that iterate through prepared earlier values and automatically changes number of output raster e.g. from outDirectRad1 to outDirectRad2.

Parameters that will change on each iteration:

e.g. first iteration

timeConfig = TimeMultipleDays(2015, 2, 32)

diffuseProp = 0.6

transmittivity = 0.2

e.g. second iteration

timeConfig = TimeMultipleDays(2015, 1, 31)

diffuseProp = 0.4

transmittivity = 0.3

**

Rest parameters are the same in each iteration

**

Here is my code:

import arcpy
from arcpy import env
from arcpy.sa import *

env.workspace = "C:/xxx/output"

arcpy.CheckOutExtension("Spatial")

inRaster = "C:/xxx/output/Test"
latitude = 47.04
skySize = 200
timeConfig = TimeMultipleDays(2015, 2, 32)
dayInterval = 14
hourInterval = 1
zFactor = 1
calcDirections = 32
zenithDivisions = 8
azimuthDivisions = 8
diffuseProp = 0.4
transmittivity = 0.6
outDirectRad = "C:/xxx/output/outDirectRad1"
outDiffuseRad = "C:/xxx/output/outDiffuseRad1"
outDirectDur = "C:/xxx/output/dir_dur1"

globRad1 = AreaSolarRadiation(inRaster, latitude, skySize, timeConfig,
dayInterval, hourInterval, "NOINTERVAL", zFactor, "FROM_DEM",
calcDirections, zenithDivisions, azimuthDivisions, "UNIFORM_SKY",
diffuseProp, transmittivity, outDirectRad, outDiffuseRad, outDirectDur)

globRad1.save("C:/xxx/output")


texts: How does Jenks optimisation method assist in ESRI mapping?

texts:
Why would someone use the Jenks optimisation method in ESRI mapping?

What would be easily located or discovered with this analysis?


texts: How to upload a NETCDF to Geoserver using cURL?

texts:
I'm trying to upload a NETCDF file (.nc) to Geoserver Store of Coverage ImageMosaic
The command is this:

curl -v -u admin:geoserver -XPOST -H "Content-type: text/plain"
-d "file:///data/coverages/myNetcdf.nc"
"http://410.152.23.128:8080:8080/geoserver/rest/workspaces/mygis/coveragestores/...
/external.imagemosaic"

as shown in Geoserver 2.4.0 documentation, but I get this message:

* About to connect() to 410.152.23.128 port 8080 (#0)
*   Trying 410.152.23.128...
* Adding handle: conn: 0x94df50
* Adding handle: send: 0
* Adding handle: recv: 0
* Curl_addHandleToPipeline: length: 1
* - Conn 0 (0x94df50) send_pipe: 1, recv_pipe: 0
* Connected to 410.152.23.128 (410.152.23.128:8080) port 8080 (#0)
* Server auth using Basic with user 'admin'
> POST /geoserver/rest/workspaces/mygis/coveragestores/.../external.imagemosaic
HTTP/1.1
> Authorization: Basic YWRtaW46TYVvc2VydmVy
> User-Agent: curl/7.32.0
> Host: 410.152.23.128:8080
> Accept: */*
> Content-type: text/plain
> Content-Length: 31
>
* upload completely sent off: 31 out of 31 bytes
< HTTP/1.1 400 Bad Request
* Server Apache-Coyote/1.1 is not blacklisted
< Server: Apache-Coyote/1.1
< Content-Type: text/plain
< Transfer-Encoding: chunked
< Date: Wed, 25 Sep 2013 15:48:37 GMT
< Connection: close
<
Failed to locate the input file file:/data/coverages/myNetcdf.nc* Closing connectio
n 0

the file EXISTS instead,
do u know how to fix this issue? is there another way to upload NetCDF files in ImageMosaic Stores? Plugin of Geoserver 2.4.0 doesn't select my file
thanx in advance


texts: Converting character attribute to numeric attribute

texts:
I am following the very simple code found in this question:

Converting shapefile to raster: Character attributes

to create a new numeric variable based on a factor. My problem is that the new numeric attribute has a complete different series of numbers than the factor column.

For example, instead of reading: 18001, 18002, 18003, it converts these numbers into 786, 894, 895: I have tried as.numeric and as.integer, and none of these work.

For further details, the shapefile column is described as:

commune_5m@data$NUMERO: Factor w/ 1842 levels "18001", "18002", ...: 786 894 895


texts: Import Map created in OCAD in ArcGIS Pro

texts:
I have a map created with the OCAD Map Builder (https://www.ocad.com/en/). I would like to import this map in ArcGIS Pro so that I can create a version of the map that is compatible with the ArcGIS Runtime for mobile Devices.

OCAD offers multiple export options including exporting to Shapefiles or DXF. Unfortunately the color and symbology information is lost when exporting in either a Shapefile or DXF.

Is it possible to somehow transfer the Symbology and color Information from the OCAD File to ArcGIS Pro?

I found the OpenOrienteering Mapper (https://github.com/OpenOrienteering/mapper) which can read OCADs .ocd files. But unfortunately it does not offer a way to export the map to any other compatible Format.

Any ideas how I could transfer as much information as possible from the OCAD file into ArcGIS Pro??


texts: QGis points distance from shore in metres

texts:
I have a grid of depth points in the Irish Sea that i'm trying to add other data fields to, in this case distance from shore. I've tried both of these suggestions to no avail, but am having some success with this one. My problem is that the resulting raster  has distances as decimal degrees rather than metres.

I've changed my project properties > canvas units to metres, re-done the rasterise, re-done the proximity, but it's still decimal degrees. Any ideas how I can get the results in metres?


texts: Calculating combined area of multiple polygons by counting overlapping area only once in R?

texts:
I'm using a shapefile in which polygons represent the geographic distribution of species. I need to calculate the combined range of several species sets, but the overlapping areas should be counted only once. My initial solution was to dissolve the polygons of the set of species of interest and then calculate its area. However, I must do this with more than 2000 sets of species and some of them have up to 10000 polygons, making the analysis take too long.

Is there another faster or more practical way of calculating the combined area of a set of polygons by counting overlapping areas only once?

For example, is there any way to delineate the outer contour of overlapping polygons?

I'm conducting these analyzes in R, but I'm willing to use other software if it decreases computational time considerably.


texts: Calculate position to put box into another container?

texts:
there has a big box, and the big box already contains many small box, then should place a new small box into the container.

How do I calculate a big enough space to place the small box?

There should be many place to hold the box.

How do I calculate at least one place?

Could not be rorate the small box.


texts: R find length along line from origin to another point

texts:
In R, I have a collection of points (~15,000) on the boundary of a polygon. For every point, I want to calculate the distance from the origin vertex (1st plotted vertex of the polygon).

I have a method that works (below) - using the lengthToPoint function defined here (with the print(paste("Found...) line commented out) - but I cannot find a quick function for distance along a line for a Spatialpointsdataframe (or whatever format) or for more than 1 point at a time.

require(maptools)
require(sp)
require(rgeos)

# create a dummy polygon: coordinates, promote to Polygons and SpatialPolygons
coordsPol = cbind(c(1,2,3,4,4),c(3,2,2,1,3))
verts <- rbind(coordsPol, coordsPol[1, ])
pol <- SpatialPolygons(list(Polygons(list(Polygon(verts)), ID="a")))

# create spatialpointsdataframe (they all fall inside the polygon)
p = data.frame(x=c(1.5,2.5,2.75,4,3.5),y=c(2.5,2,2,1.3,3),ID=c("a","b","c","d","e"))
coordinates(p) <- ~x+y

# the polygon needs to be a SpatialLines object for the lengthToPoint fucntion
l <- as(pol,"SpatialLines")

# the SPDF needs to be SpatialPoints for the lengthToPoint function
pm <- as(p,"SpatialPoints")

# create a vector of distances from the origin using lengthToPoint as this can
# only take one point at a time
# https://gis.stackexchange.com/questions/209254/calculate-distance-of-points-spatialpoints-object-along-a-path-spatialline
lenvec <- vector()
for(pp in 1:length(pm)){
  lenvec[pp] <- lengthToPoint(l,pm[pp])
  }

p@data$dist <- lenvec

Is there a function that returns the distance along a line from the origin for > 1 point at a time?


texts: Is there any fast looping through too many shapefiles in a folder using ArcPy?

texts:
Looping through shapefile using arcpy.ListFeatureClasses() is very slow for many shapefiles. because first thing first, it writes all in memory as a list and then read each shapefile one by one.
Is there any other way to read one by one shapefile without making them as a list and write it in the memory:
This is the slow script:
arcpy.env.workspace = shapefile_folder
fc_list = arcpy.ListFeatureClasses()
for shapefile in fc_list:
     print (fc) # doing the processing over shp here


texts:
I am trying to run a species distribution model using random forests in R.
I have location data for species occurrence in a national park.
When I load the WorlClim data, the data for my coordinates is all NA.

Why is this?

Is it because my overall location is too small/narrow?

When I ask for random WorldClim data in my location is gives it?

> r <- getData("worldclim",var="bio",res=10)
> r <- r[[c(1,12)]]
> values <- extract(r,points)
> df <- cbind.data.frame(coordinates(points),values)
> head(df)
  lon       lat Temp Prec
1 xx.16025 -xx.42106   NA   NA
2 xx.15480 -xx.41602   NA   NA
3 xx.15359 -xx.18017   NA   NA
4 xx.14291 -xx.26207   NA   NA
5 xx.13306 -xx.14253   NA   NA
6 xx.12042 -xx.00593   NA   NA

> e <- extent(SpatialPoints(data))
> bg <- sampleRandom(r, 5000, ext=e)
> head(bg)
     Temp Prec
[1,]  208  460
[2,]  203  482
[3,]  198  491
[4,]  201  472
[5,]  199  473
[6,]  194  490


texts:
I'm trying to create a map for a fictional map. I'm using Leaflet to display the map, with tiles, different levels of zoom and so on.

Now I need to generate, to "draw", some vector data. Specifically, I have streets that I need to add to the map.

Is there any software where I can do something like "draw" those streets as a vector and export them to use in a JS file?


texts: Openlayers2 Vector Feature label not Visible in IE 8 but visible in Google Chrome....?

texts:
I have created a WEB MAP using Openlayers-GeoServer-PostGIS....

The Label of one of the Vector feature layer is not visible in Internet Explorer 8(IE 8) ,but clearly visible in Google Chrome 27.0.1453.94....

As most of my map viewers use IE8 so it is important to have the full map along with label in IE8...

part of the code used is as follows...

var vectorLayer = new OpenLayers.Layer.Vector("SupplyOff");
    map.addLayer(vectorLayer);

<?php

$qry_xy="select * from ci_master@ltdvlp  where  x is not null  AND GIS_SOURCE_CODE LIKE '01%' AND ROWNUM < 5";

$parsed =oci_parse($db_conn_prod,$qry_xy);
oci_execute($parsed);
$nrows3=oci_fetch_all($parsed,$results3);
echo $nrows3;

for($i=0;$i<$nrows3;$i++)    {
$x=$results3["X"][$i];
$y=$results3["Y"][$i];
$consno=$results3["CONS_NUM"][$i];
$source=$results3["GIS_SOURCE_CODE"][$i];
  ?>

var description_new = '<?php echo $consno; ?>';
            var x = '<?php echo $x;  ?>';
            var y ='<?php echo $y;  ?>';
        var cons ='<?php echo $consno; ?>';
            var type = '<?php echo $source; ?>';

    vectorLayer.addFeatures([new OpenLayers.Feature.Vector(
                new OpenLayers.Geometry.Point(x,y ),

                                  {
                    //Attributes go here
                    description:'cons: '+cons,
                    size: 5 + (Math.floor(Math.random() * 20)),
                    label: 'F'+cons,
                    strokeWidth: (Math.floor(Math.random() * 10))
                }
            )]);

<?php

            }
?>

    //Create a style object to be used by a StyleMap object
    var vector_style = new OpenLayers.Style({
        'cursor': 'pointer',
        'fillColor': '#2bff05',
        'fillOpacity': .8,
        'fontColor': '#a52505',
        'label': '${label}',
        'pointRadius': '${size}',
        'strokeColor': '#232323',
        'strokeDashstyle': 'solid',
        'strokeWidth': '${strokeWidth}',

    });

    //  style object we just created
    var vector_style_map = new OpenLayers.StyleMap({
        'default': vector_style
    });

    //Add the style map to the vector layer
    vectorLayer.styleMap = vector_style_map;

map.addLayer(vectorLayer);
var controls = {
      selector: new OpenLayers.Control.SelectFeature(vectorLayer, { onSelect: createPopup, onUnselect: destroyPopup })
    };

    function createPopup(feature) {
      feature.popup = new OpenLayers.Popup.FramedCloud("pop",
          feature.geometry.getBounds().getCenterLonLat(),
          null,
          '
'+feature.attributes.description+'
',
          null,
          true,
          function() { controls['selector'].unselectAll(); }
      );

      map.addPopup(feature.popup);
    }

    function destroyPopup(feature) {
      feature.popup.destroy();
      feature.popup = null;
    }

    map.addControl(controls['selector']);
    controls['selector'].activate();


texts:
I have a geotiff and I have the corner coordinates set using gdal_translate and the a_ullr argument setting out the min/max long/lat. This provides a bounding box for the geotiff to be projected.

But I want this projected different to just a bounding box when displayed as a layer on a map. E.g. see this picture below.

On the left is what the geotiff is, with the corner coordinates defined. But I want each corner to have it's own lat/long defined separately, so you get as shown on the right.

Whats the best way of doing this so it's displayed as so as a layer on a map.

(Using ESPG:4326)


texts: How do you change the extent of a clipped raster layer in qgis to be the mbb?

texts:
I have a raster layer that I have clipped using a polygon in qgis. The extent of the raster layer did not change.
How do I update the extent of my raster layer so that it is the minimum bounding box of my clipping?


texts: How to sort field by ascending numerical order in QGIS 2.14?

texts:
I work with this field (in categorized style):

I try to find how can I sort the "value" column (in the style properties) by ascending numerical order. In the attribute table "value" is "YEUD" field:

"YEUD" field id numerical field:

The problem is when I click on "value" name, it sorts the field not the way I am seeking.

I read Why are numeric values in attribute table not sorted correctly? and my version is earlier,  and also I don't want to update my QGIS version.


texts: How to compile QGIS 2.4.0 on Ubuntu Precise 12.04

texts:
Been trying to compile this for a couple of hours, using this download:

http://qgis.org/downloads/qgis-latest.tar.bz2

and these instructions:

http://htmlpreview.github.io/?https://raw.github.com/qgis/QGIS/master/doc/INSTALL.html#toc3.

I have problems with libproj-dev - I can't install it, get this error:

  The following packages have unmet dependencies:
  libproj-dev : Depends: libproj0 (= 4.8.0-3~precise2) but 4.8.0-3~precise5 is to be installed
  E: Unable to correct problems, you have held broken packages.

So I installed the headers by hand - compiler accepted them - and have tried to link in /usr/lib/libproj.so as PROJ_LIBRARY in ccmake. The compiler doesn't like this, I get:

  make[2]: * No rule to make target /usr/lib/libproj.so', needed byoutput/lib/libqgis_core.so.2.4.0'.  Stop.
  make[1]: * [src/core/CMakeFiles/qgis_core.dir/all] Error 2

I also get 'Could not find GRASS' a little earlier in the build - don't know if this is significant.

I'm a bit confused, libproj.so shouldn't be a target? Anyway, can anyone help? My ppa is the ubuntugis stable one - maybe I need unstable?


texts:
I am using QGIS 3.4 LTS on Kubuntu Bionic LTS. I am using several global datasets like NaturalEarth, HydroRIVERS etc. Loading all of them takes too much time and memory especially on my laptop.

I am mostly working in India area only. So I'd like to apply some settings in my project to tell QGIS to load features from within India only for performance reasons. It is acceptable to have it as a lat/lon extent constraint rather than clipping to a polygon.

There is existing advice on how to clip data to given polygons but I don't want to maintain separate shapefiles etc with clipped data. So I'd prefer to set it in the project using global upstream data.

How can I achieve this in QGIS?

This question may be the QGIS equivalent of Definition query to display only part of layer that overlaps with specific polygon using ArcGIS Desktop.


texts: Automatically set PostGIS/QGIS style to a layer

texts:
I have a plugin for my users which allow them to create a spatial view for different layers. I have define a set of styles fitting the display I want to give. I would like to define these styles to the view automatically every-time a new view is created. I am currently using QGIS 2.14, PostGIS 2.3, PostgreSQL 9.5 and Python 3 with the psycopg2 library.

I thought about simply populating the layer_style table in the database, but without success.

Does anyone know how to do so ?


texts: How to use cartodb.js without storing sql in javascript

texts:
Using cartodb.js, how can I change the layer definition without building the sql dynamically in javascript and passing it as a parameter.  Aside from revealing my table structure (not a huge deal at the moment, but it might be later on), it seems like a bit of a maintenance headache.

mylayer.set({
  sql: "select * from mytable WHERE somefield = 'x'"
});

I don't see (but could be missing) options within CartoDb itself.  Another option is I could set up a separate web server where I would send a request (e.g. someotherserver.com/request?somefield=x) which could then generate sql and then talk to CartoDB and return the results.  The problem with that approach is that calling layer.set(...) provides a good deal of additional functionality that I'm not sure my separate web server would be able to provide.

Does anyone have a good solution to this problem?


texts: Merging boundaries with boundary shapefile in ArcMap?

texts:
I have a boundary shapefile of a city with 19 district boundaries, and I need to reflect some changes and create 16 districts instead.

Boundaries I want to merge are all neighboring districts (two neighboring districts need to become one).

I'm a beginner so I'm not sure what I should be doing to get at this, please give me some advice! I'm using ArcMap 10.5


texts: How to make GRASS 7 addon available in QGIS 2.12?

texts:
I stumbled across the GRASS-Tool i.segment but could not find it in my QGIS processing toolbox. I found out that it is an addon for  GRASS 7 that I already installed in GRASS 7 standalone. How can I import Grass 7 into QGIS?

I installed the latest QGIS Version with the OSGeo4W installer and enabled Grass for installation. The thing is, that the OSGeo4W installation installed Grass 6.4.3 and not Grass 7.  I do know how to configure external applications (I have already done this for Grass 6.4.3 and for OTB). I am wondering how it is possible to make GRASS 7 (tools and addons) available within the processing toolbox?

I updated Windows 7 to Windows 10.


texts: ArcGIS labels with expression not appearing

texts: I have a layer which should be labelled based on two different fields. I've chang the Expression in Layer Properties from [FIELD1] to [FIELD1] + '\n' + [FIELD2]. However, only Field1 is showing in my map, never Field2. I have fixed the scale range and I am using the Maplex Label Engine. Why aren't both fields showing in my label?


texts: Are there any tools which can autopopulate multiple features with metadata fields in ArcGIS?

texts: I have an enterprise geodatabase which contains hundreds of feature classes.  I would like to publish the metadata from each feature to geoportal server in FGDC format.  Is there an out of the box tool/script which allows users to pre-populate fields for all data sets in a geodatabase?  Ideally, features such as originator, date, etc.


texts: How to make focus on Leaflet marker with a click from another window

texts:
I am making a project using Leaflet. I made a script which creates a map with one marker. A click on the map creates a window with the marker coordinates and I would like to be able to click on coordinates in the new window and having a focus on the marker linked to this coordinates on the map. But I have no idea on how to do it or if there is a better way to do this even after searching on Internet. You can easily launch the code.

<!DOCTYPE html>

    <title>Leaflet</title>

    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
     <link rel="stylesheet" href="https://unpkg.com/leaflet@1.5.1/dist/leaflet.css"
   integrity="sha512-xwE/Az9zrjBIphAcBb3F6JVqxf46+CDLwfLMHloNu6KEQCAWi6HcDUbeOfBIptF7tcCzusKFjFw2yuvEpDL9wQ=="
   crossorigin=""/>

<div


texts:
I have:

GeoServer
Oracle Spatial
WMS

To communicate with Oracle, GeoServer uses the schema, username and password.
Is it possible to use multi-user to communicate with Oracle locator table or switch the Oracle?


texts: I am very new to QGIS and I simply have two tif images which I need to do pixel analysis on so I can compare the change in vegetation over time. I have no idea how to go about this and if anyone has some advice for me I would be very grateful


texts: Why maximum value is always same in pyqgis?

texts:
I want to get the maximum value of the layer. (ex, 1st field)

I used the layer (or provider) maximumValue () function.

but the same value is always obtained.

ex>

There are five features and the first field value is 1, 2, 3, 4, 5.

The maximumValue () function returns 5.
Delete the fifth feature. The remaining features are four, and the first field value is 1,2,3,4.
The maximumValue () function returns 5 equally.

I expected it to return 4.

Why does the same value get returned after the feature changes?

(commitChanges () was done.)

aLayer = QgsMapLayerRegistry.instance().mapLayersByName('aaaLayer')[0]
print aLayer.getValues('field_1')
>>> ([1, 2, 3, 4, 5], True)  --------------> field_1 value list : 1, 2, 3, 4, 5
print aLayer.maximumValue(0)
>>> 5 -----------------------------> 1
aLayer.startEditing()
aLayer.deleteFeature(5) # ---------------> delete feature.
aLayer.commitChanges()
print aLayer.getValues('field_1')
>>> ([1, 2, 3, 4], True)  ------------------> field_1 value list : 1, 2, 3, 4
print aLayer.maximumValue(0)
>>> 5  ------------------------------> 2

Why are 1 and 2 the same?

How can I solve it?


texts: I would like to export from QGIS to another map (eg umap, openstreetmap) in kml or other format of the graduated points with different measures based on the concentration of certain events. Unfortunately, exporting them to any format I find only the markers on the destination map. Do you know how I can solve?


texts: Crop a raster using rasterio and geopandas

texts:
I am cropping a set of historic aerial photographs. These photographs have large black areas in the edges (0 value). However there is also valid data with a 0 value. The workflow I am using is:

Load raster with rasterio
Polygonize the raster with rasterio.features.shapes()
Identify polygons where value = 0 and size > 5000 square meters
Mask the original images with polygons, perform an inverted mask

Here is my current code for masking a single image:

import rasterio
from rasterio import features
from rasterio import mask
import json
import geopandas as gpd

results = []
final_results = []

with rasterio.open(r"C:\1927_oahu\tif\_Line1_6to8_0.tif") as src:
    src_meta = src.meta
    src_affine = src_meta.get("transform")

    band = src.read(1)

    for geometry, raster_value in features.shapes(band, transform=src_affine):
        if raster_value == 0:
        result = {'properties': {'raster_value': raster_value}, 'geometry': geometry}
        results.append(result)

gpd_results = gpd.GeoDataFrame.from_features(results)

gpd_results["area"] = gpd_results["geometry"].area

gpd_results_filtered = gpd_results[gpd_results["area"] > 5000]

gpd_filtered_json_str = gpd_results_filtered.to_json()

gpd_filtered_json_dict = json.loads(gpd_filtered_json_str)

for k, v in gpd_filtered_json_dict.iteritems():
    if k == "features":
        for items in v:
            #final_results = {"coordinates": (items.get("geometry").get("coordinates"))}
            final_results = {"geometry": (items.get("geometry").get("coordinates"))}

masked_band, masked_transform = mask.mask(src, final_results, invert=True)

src_meta.update(dtype=rasterio.uint8, nodata=0)
with rasterio.open(os.path.join(r"C:\1927_oahu_output", "out.tif"), 'w', **src_meta) as dst:
    dst.write_band(1, masked_band.astype(rasterio.uint8))

When I run this code I receive the following error: AttributeError: 'str' object has no attribute 'get'

The documentation for rasterio.mask states: Polygons are GeoJSON-like dicts specifying the boundaries of features in the raster to be kept. All data outside of specified polygons will be set to nodata.

I am assuming that I am giving rasterio.mask the wrong type of "GeoJSON-like dict". I have tried to reformat the dict several ways without success. Does anyone know the correct way to convert GeoJSON to a "GeoJSON-like dict"?

Or can anyone provide the correct format of a "GeoJSON-like dict"?

I am new to rasterio and geopandas.


texts:
I have a selection in vector point. I would like to do random sample from it. I am using qgis to do it.

Here is my selection code:

selection=layers_canvas[0].getFeatures(QgsFeatureRequest().setFilterExpression(u'"fitofissio" ILIKE \'%Vereda%\''))
layers_canvas[0].setSelectedFeatures([s.id() for s in selection])

I don't know how to do it. The result of select is on class and it isn't able to use ramdom.sample in it. I tried to do my random sample with this function.


texts: Advanced label expression in ArcMap to exclude part of a string

texts:
How do I write an advanced labeling expressions to only show a part of a string?

I want my data that looks like this EVA007 only to show the 007 part.


texts: I'm interested in accessing a shapefile or tabular dataset with household income data at as granular of level as possible in Canada (subunits on the island of Montreal specifically). Is this data available at all in Canada? Or would there be a similar dataset that would speak to socioeconomic levels at fairly fine scales (e.g., town or parts of a town)?


texts: How can I programmatically project the Z value of a z-aware feature point, using ArcObjects in .NET/C#?

texts:
I discovered the following code leaves the Z coordinate of featurePt untouched:

// if not WGS84 project it
if (!SpatialRefsEqual(featureSpatialRef, srWgs84))
{
    // featurePt here is a Z-Aware instance of IPoint
    // use the Geometry Interface to do the actual project
    featurePt.SpatialReference = featureSpatialRef;
    featurePt.Project(srWgs84);
}

I see the IVerticalCoordinateSystem interface but it unclear to me how to use it to do a projection.  For example, simply plugging it in via:

featurePt.Project(myInstanceOfIVerticalCoordinateSystem)

still leaves Z untouched.

It also seems like I could just write a routine to convert the Z-coordinate in the linear units of the layer to my desired linear units (meters in this case).

Is there a best practice for projecting the Z-coordinate?

EDIT

According to mkennedy at ESRI, ArcGIS does not provide a true vertical transformation (thank you for answering!).  If anyone knows of any outside library to vertically transform among vertical datums please let me know.


texts: Is it possible to change number format to scientific in QGIS layer legend?

texts:
I have a layer with really small values and the classes in the legend do not look nice both for work in QGIS but also in the print layout.

Is it possible to somehow change the format of the numbers in the layer legend so instead numbers like:
0.000012345
I could get scientific notation:
1.2345e-5
or do I have to manually convert the numbers from the legend somewhere "outside" - maybe in some online tool like this and then paste it back in the legend as class label?
I googled and also looked in the QGIS documentation but did not find any solution.


texts: Annotation/Label size adjusting to map scale in QGIS

texts:
I'm trying to label the street names for a map I'm working on but I do not want the labels to automatically adjust to the map scale when zooming in or out. My question is: In QGIS how do I set a fixed scale for the labels so that it wont change with the map scale?
In the labeling properties I have checked the box "Scale-based visibility" under the rendering tab but it doesn't seem to affect it at all.


texts: Automatically check geographic coordinates?

texts:
I have a bunch of localities worldwide for which I need to check their geographic coordinates. So far, I only count on information on the corresponding countries, localities' names, and "suggested" latitudes and longitudes in decimal degrees. I'm wondering if is there any 100% free tool (ideally an R package) that allows me to automatically check the degree of error of some coordinates.
Ideally, I'm thinking about a system that (1) picks the information available, (2) checks it in some gazetteer or similar the localities around the provided coordinates (let's say, within a buffer of 20 km) and returns if there is a locality with a similar name within the buffer.
Any idea?


texts: Calculating summary statistics and access resulting value in table using ArcPy?

texts:
With arcpy, I want to calculate the mean value of a specific field in a shapefile and print the resulting value for the user.
I use the Summary Statistics tool and it works, but the resulting value is stored in a table at a specific path, but how can I access the calculated value with python?

Here is my code. It will only output the path to the table, but I need the calculated value.

stats = arcpy.Statistics_analysis(myshapefile, "G:/temp/stats", [[field, "MEAN"]])
arcpy.AddMessage(stats.getOutput(0))


texts: Where could I download Income data at the Block Group level... census data

texts: Where could I download  Income data at the Block Group level,   actually I'm not sure where or how to download it at the tract level either.    I would be interested in the most recent data that is available.


texts: Where can i find Latitude & Longitude data for the coastlines of the world?

texts:
Im developing an application to print flight paths on a world map.

To do this i need a dataset to describe the latitude and longitude of various coastlines (and country borders if possible) of the world so I can plot them.

Is there anywhere on the web I can find this kind of data? (preferably free/cheap)

Thanks!


texts: Determining ArcCatalog's currrently selected folder using ArcPy

texts:
Is there a way to get the current directory opened in ArcCatalog 10.0 using ArcPy?

Something like os.getcwd() but instead gives the selected folder on the catalog tree.


texts: Using "Crtl + C" to copy an item from QTableWidget in QGIS, but get all attributes from a feature

texts:
I am using PyQGIS to develop a table as a QGIS 3.4 plugin.
The table is created by QTableWidget and looks good.
The weird thing is when I try to copy (Ctrl+C keyboard shortcut) an item in the table, it indeed copies the whole feature, like Copy Feature action in the Identify Features tool.
Basically, the item is just a simple string. I really don't know why this behaviour happened.
e.g. in the screenshot below, I only copied the blue highlighted item. But when pasting, I got an output of all attributes of the corresponding feature.
I didn't implement any function related to the Copy functionality.

Same result as I can get from Copy Feature action.

My Script:
    self.ui.scrollArea.setWidgetResizable(True)
    self.ui.widget = QtWidgets.QWidget()
    self.ui.scrollArea.setWidget(self.ui.widget)
    self.ui.layout_VArea = QtWidgets.QVBoxLayout(self.ui.widget)
    font = QtGui.QFont()
    font.setBold(True)
    font.setItalic(True)
    font.setPointSize(9)

    # use QtableWidget to create table
    self.ui.tableWidget = QtWidgets.QTableWidget(self.ui.widget)
    self.ui.tableWidget.setObjectName("tableWidget")
    self.ui.layout_VArea.addWidget(self.ui.tableWidget)
    self.ui.tableWidget.setColumnCount(len(self.userFields[self.layer]))

    labels = [self.labels[self.layer].get(field, field) for field in self.userFields[self.layer]]
    print(labels) # a list of string

    self.ui.tableWidget.setHorizontalHeaderLabels(labels)
    print(self.ui.tableWidget.selectionBehavior()) # 0
    print(self.ui.tableWidget.selectionMode()) # 3
    # get how many rows in the values
    rows = len(values)

    # set table rows
    self.ui.tableWidget.setRowCount(rows)

    for indexRow, id in enumerate(values.keys()):
        rowString = list(map(str, values[id]))
        print(rowString) # a string list
        for indexColumn, item in enumerate(rowString):
            flags = Qt.ItemIsSelectable | Qt.ItemIsEnabled
            qItem = QtWidgets.QTableWidgetItem(item)
            qItem.setFlags(flags)
            self.ui.tableWidget.setItem(indexRow, indexColumn, qItem)


texts: Why [0] is needed in row[0] when using ArcPy (Data Access) cursors?

texts:
when running the following 2 sample codes, "[0]" is needed to be added after "row" in order to make code run well. I was confused. can anyone help explain? I know that it might be too easy for most of you from here.

sample code #1:

sf = r"C:\Data\Exercise07\alaska.shp"
cursor = arcpy.da.UpdateCursor(sf,["PERIMETER"])
for row in cursor:
...     if row[0] == 0.224:
...         print "find"

sample code #2:

cursor =arcpy.da.SearchCursor("alaska",["PERIMETER"],'"PERIMETER"=0.224')
for row in cursor:
...     print row[0]


texts:
I want to add a classification to photos which are in an excel spreadsheet within a KML file.

Need a way to combine/join fields that are in an Excel spreadsheet to a KML file.


texts: How to update a value in an ArcSde table using ArcMap's VBA

texts:
First of all, I want to update a table, not a feature class (no geographical information).
I tried the following

Public Sub updateInSql(tableName As String, whereClause As String, field As String, value As String)

    ' get table
    Dim pTable As ITable: Set pTable = getTableFromSQL(tableName) ' this is my function that gets the table object (it works for sure)
    If (pTable Is Nothing) Then Exit Sub

    Dim pQueryFilter As IQueryFilter: Set pQueryFilter = New QueryFilter
    pQueryFilter.whereClause = whereClause
    If (whereClause = "") Then Set pQueryFilter = Nothing

    Dim pCursor As ICursor: Set pCursor = pTable.Update(pQueryFilter, False)
    Dim pRow As IRow: Set pRow = pCursor.NextRow
    Do Until pRow Is Nothing
        pRow.value(pRow.fields.FindField(field)) = value
        pRow.Store
        Set pRow = pCursor.NextRow
    Loop

End Sub

I get "Automation Error" at "Set pCursor = pTable.Update(pQueryFilter, False)"
Why?


texts: Getting coordinates falling inside two known points using javascript

texts:
I'm working on a solution, where a parent track the vehicle in which his kids are travelling. Initially I will be getting all the geo fence data from the system (Like stopings) and store it in an array, and I'll also be getting the current position of the vehicle every 1 minute.

Now I want to find out which is the next stop or geofence approaching in any given time, considering current position as origin and the parent location as destination, need to find out which is the next stop inbetween those two points?

These are the ideas on which I worked on

Calculate the distance from the current position to all the geo fence data and select the short distance. But it did not work out because the next stop can be far from the previously crossed geofences, so that I only get the previous geofence not the next
Using 2 points I calculated the heading angle, but with this I couldn't find the next stop
Finally I stubbled upon a usable solution, which is to use intersection based on 2 positions, calculating the geofences between lat (min), lat(max) and lng(min) and lng(max). Then finding the distance from the current point and sorting to find the next stop.

But this solution failed because, when the vehicle route is Zigzag the fence data can be outside the bounder of origin & destination points

Below is the images on this solution

Please help me in finding a solution...


texts: Google Earth Engine - Blend Sentinel-1 SAR Images?

texts:
Does anybody know why I am getting these strips across my Sentinel-1 images and how to fix it?

I am using multiple images over Mozambique between January 1st and January 13th.
I want to use this as a "before" image to then calculate flooding extend after a cyclone.
I assume I need to merge the images together or something to get median values between the 3 overlapping images?

This is the code I am using (including code to generate flooding extent in each month from March through to October.

// Filter the collection for the VV product from the descending track
var imageCollection = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterBounds(Area)
    .filter(ee.Filter.date('2019-01-01', '2019-10-20'))
    .select(['VV']);
print(imageCollection);

//create map layer of before flooding
var beforeJanuary = imageCollection.filterDate('2019-01-01', '2019-01-13');

//create map layer for every month following cyclone landfall
var afterMarch = imageCollection.filterDate('2019-03-14', '2019-03-29');
var afterApril = imageCollection.filterDate('2019-04-01', '2019-04-30');
var afterMay = imageCollection.filterDate('2019-05-01', '2019-05-31');
var afterJune = imageCollection.filterDate('2019-06-01', '2019-06-30');
var afterJuly = imageCollection.filterDate('2019-07-01', '2019-07-31');
var afterAugust = imageCollection.filterDate('2019-08-01', '2019-08-31');
var afterSeptember = imageCollection.filterDate('2019-09-01', '2019-09-30');
var afterOctober = imageCollection.filterDate('2019-10-01', '2019-10-20');

//Obtain an average pixle value for overlapping vv images
var VVJanuary = beforeJanuary.median();
var VVMarch = afterMarch.median();
var VVApril = afterApril.median();
var VVMay = afterMay.median();
var VVJune = afterJune.median();
var VVJuly = afterJuly.median();
var VVAugust = afterAugust.median();
var VVSeptember = afterSeptember.median();
var VVOctober = afterOctober.median();

//add layers to map
Map.addLayer(VVJanuary, {min: -14, max: -7}, 'Before');
Map.addLayer(VVMarch, {min: -14, max: -7}, 'March');
Map.addLayer(VVApril, {min: -14, max: -7}, 'April');
Map.addLayer(VVMay, {min: -14, max: -7}, 'May');
Map.addLayer(VVJune, {min: -14, max: -7}, 'June');
Map.addLayer(VVJuly, {min: -14, max: -7}, 'July');
Map.addLayer(VVAugust, {min: -14, max: -7}, 'August');
Map.addLayer(VVSeptember, {min: -14, max: -7}, 'September');
Map.addLayer(VVOctober, {min: -14, max: -7}, 'October');
Map.addLayer(Area)

//Centre the map view over our region of interest
Map.centerObject(roi, 20);


texts: Creating tool from script using ArcPy gives error 000714

texts:
The script is fully functional, but I need to convert it into a tool which will allow the user to customize work environment, expression, name of shapefile, name title text, update date text automatically, and select a output file.

The error message is:

  The Microsoft Access database engine could not find the object 'GDB_ReleaseInfo'. Make sure the object exists and that you spell its name and the path name correctly.

I am new to ArcPy.

   # Import arcpy module
... import arcpy
...
... #Set Work Environment Default File
... myWS = arcpy.GetParameterAsText(0)
... if myWS == '#' or not myWS:
...     myWS = r"C:/Tmp/GIS_Tmp" # Default value if none passed to script
...
... #arcpy.env.workspace = myWS
... #arcpy.env.workspace = "C:\Tmp\GIS_Tmp"
...
... # Script arguments
... Parcels = arcpy.GetParameterAsText(1)
...
... Selection_type = arcpy.GetParameterAsText(2)
... if Selection_type == '#' or not Selection_type:
...     Selection_type = "NEW_SELECTION" # provide a default value if unspecified
...
... Distance__Buff1 = arcpy.GetParameterAsText(3)
... if Distance__Buff1 == '#' or not Distance__Buff1:
...     Distance__Buff1 = "300 Feet" # provide a default value if unspecified
...
... Distance__Buff2 = arcpy.GetParameterAsText(4)
... if Distance__Buff2 == '#' or not Distance__Buff2:
...     Distance__Buff2 = "500 Feet" # provide a default value if unspecified
...
... Expression = arcpy.GetParameterAsText(5)
... if Expression == '#' or not Expression:
...     Expression = "\"AIN\" = ''" # provide a default value if unspecified
...
... v300ft_shp = arcpy.GetParameterAsText(6)
...
... v500ft_shp = arcpy.GetParameterAsText(7)
...
...
... # Local variables:
... #Parcel_Address = Parcels
...
... #If layer already exist overwrite
... arcpy.env.overwriteOutput=True
...
... #Create Parcels layer
... arcpy.MakeFeatureLayer_management("C:/Tmp/GIS_Tmp/Parcels.shp","Parcels.lyr")
...
...
... # Process: Select Layer By Attribute
... arcpy.SelectLayerByAttribute_management("Parcels.lyr", Selection_type, Expression)
...
... # identify the current map document
...
... mxd = arcpy.mapping.MapDocument('CURRENT')
...
... # identify the data frame
... df = arcpy.mapping.ListDataFrames(mxd)[0]
...
... # use the zoom to selected features method of the data frame to update the extent
... df.zoomToSelectedFeatures ()
...
... #Set Scale
... df.scale = 3860
...
... # Process: Buffer
... arcpy.Buffer_analysis("Parcels.lyr", "C:/Tmp/GIS_Tmp/300ft_shp", Distance__Buff1, "FULL", "ROUND", "NONE", "", "PLANAR")
...
... #Import Symbology layer to buffer
... arcpy.ApplySymbologyFromLayer_management ("300ft_shp", "C:/Tmp/GIS_Tmp/Symbol_300ft.lyr")
...
... # Process: Buffer (2)
... arcpy.Buffer_analysis("Parcels.lyr", "C:/Tmp/GIS_Tmp/500ft_shp", Distance__Buff2, "FULL", "ROUND", "NONE", "", "PLANAR")
...
... #Import Symbology layer to buffer
... arcpy.ApplySymbologyFromLayer_management ("500ft_shp", "C:/Tmp/GIS_Tmp/Symbol_300ft.lyr")
...
... # Update Text within Title
... import arcpy
... mxd = arcpy.mapping.MapDocument("CURRENT")
... for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
...     # Text not within Element
...     if elm.text == "Text":
...         elm.text = "fine"
... arcpy.RefreshTOC()
... del mxd
...
... #Update Text for Date
...
... from datetime import date, timedelta
... mxd = arcpy.mapping.MapDocument("CURRENT")
... #Where my code says "DateTimeElm" your text element needs to be named the same or change this name to reflect your text element
... elm = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "DateTimeElm")[0]
... elm.text = date.strftime(date.today() - timedelta(days=1), "Current: " + '%B %d, %Y')
... arcpy.RefreshActiveView()


texts: How can I optimize ArcMap Editor's streaming mode for use with a tablet such as the Wacom Intuos?

texts:
I have the Wacom Intuos 4 Large tablet which I use occasionally outside of ArcGIS. I now need to digitize a large volume of polygons from aerial imagery and am finding the tablet's ArcGIS functionality cumbersome.

The pen tip works fine in a regular mode where I tap it to create a vertex (albeit far too easy to mistakenly engage a double click) but I would like to be able to delineate features as if I were doing it on paper - a fluid hand motion.  I turned on streaming and set it to 5 map units which is OK but when I reach a point where I need to lift and reposition my hand on the tablet or pan the page while still digitizing, streaming does not stop when I lift the pen from the tablet.  Ideally, I would like the streaming to behave exactly as it does in Google Earth - streaming pauses when the tip is lifted and resumes when a pressure is again applied to the tip. GE also generalizes sketches on completion not while sketching as ArcMap does - this makes the finished product simple but precise.

Alternatively, I used ArcMap's Tablet toolbar which streams beautifully to produce "ink" sketches. However, when I lift the pen from the surface, line is interrupted and while this is fine for polylines, the line interruptions create artifacts when creating polygons.

How to I set up ArcMap Editor streaming so that streaming pauses when I lift the pen tip from the tablet's surface?

OR

How can I draw an ink sketch and pause drawing when the pen is lifted and have the sketch resume in a new location where the pen touches the surface?

OR

How can I draw an ink sketch and snap to another ink sketch start/end node?

I am off to give it this a try in QGIS...


texts: How to mark a parallel path of existing line in JOSM/QGIS?

texts:
I have a railway line imported from OSM to JOSM. This is converted as map tile (using maperitive) for web mapping in leaflet. Now I want to mark a parallel line to the existing railway line.
For example. A railway line (curved) is having a distance of 10 km. Say from 4th km to 6th km if some work has been done, it should be marked with different colour in the web map.
How to do it?


texts:
I am using spatialite and it does not contains the method ST_dumppoints and ST_dump. I want to divide several linestring into smaller parts. I also have points that I combined and that resulted in a linestrings. Now I want to divide these linestring into several parts.

I have tried using ST_line_substring and ST_line_locate_point but nothing is working.

Is there a way to divide linestring into smaller segments say 100m in spatialite?

What I have done is:

select
ASText(ST_Line_Substring(linearea.geometry,0,
ST_Line_Locate_Point(linearea.geometry, Transform(setsrid(pointarea.geometry,4326),3857))))  AS Geometry
from LineArea linearea,PointArea pointarea

PontArea is Area which I combined to make LineArea using MakeLine(ST_point(long,lat))


texts: QGIS - Installation issue Windows 7 64 bit

texts: I am currently trying to install QGIS but it has been on the same point of the installation for almost 30 minutes. I don't know if it is simply taking a long time to create a folder, but it seems as though there is an error in the installation process. It appears as if it's in the "finishing" process, but there's a bug.


texts:
I'm working through some examples in the Python Geospatial Development book by Erik Westra.

On page 439, about setting up a tile map server, he gives the following code for calculating a bounding box based on a centre location and a zoom level:

TILE_WIDTH = 256
TILE_HEIGHT = 256

x_extent = _unitsPerPixel(zoomLevel) * TILE_WIDTH
y_extent = _unitsPerPixel(zoomLevel) * TILE_HEIGHT

min_long = x * x_extent - 180.0
min_lat = y * y_extent - 90.0
max_long = min_long + x_extent
max_lat = min_lat + y_extent

def _unitsPerPixel(zoomLevel):
    return 0.703125/math.pow(2, zoomLevel)

Nowhere does he specify what the x and y units are expected to be. I assumed lat long coords (same as map) but then the above calculations do not work.

For example, suppose I have x = 50.0 and y=-2.09 and a zoom level of 5.

For the min_lat, I will then have approximately:

min_lat = -2.09 * 5.63 - 90.0

Which is going to give less than -90, the limit of our map...?!?

Everytime you pass a negative lat/long coordinate this will fall down....

Just what x/y coordinates are being expected here?


texts: How can I get the path for an input layer in a QGIS 3 script?

texts:
Writing scripts in previous versions of QGIS, I would get the layer path to pass to algorithms using processing.getObject(layer). In QGIS 3, the scripting system has been rewritten, and I now define a 'source' for input layers like this:

SPECIES_POINTS = 'SPECIES_POINTS'      # INPUT = 'INPUT'

[...]

speciesPointSource = self.parameterAsSource(
            parameters,
            self.SPECIES_POINTS,
            context
        )

I then try to pass it to an algorithm:

pointsReprojected = processing.run("native:reprojectlayer", {
            'INPUT': speciesPointSource,
            'TARGET_CRS': CRS,
            'OUTPUT' : speciesExtent_id
        }, context=context, feedback=feedback)['OUTPUT']

but when I try to run this, I get an error:

Unable to execute algorithm
Could not load source layer for INPUT: invalid value
Traceback (most recent call last):
File "<string>", line 191, in processAlgorithm
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\tools\general.py", line 96, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\core\Processing.py", line 142, in runAlgorithm
raise QgsProcessingException(msg)
_core.QgsProcessingException: Unable to execute algorithm
Could not load source layer for INPUT: invalid value

I can see that speciesPointSource has type QgsFeatureSource, so I'm not surprised this doesn't work, but I can't see any attributes that would give me the file path for the layer. I have also tried self.SPECIES_POINTS but then I get the error Could not load source layer for INPUT: SPECIES_POINTS not found. What should I be passing as 'INPUT'?


texts:
I have a polygon and known 2 points. One is the center of the polygon and another is outside the polygon. I'm trying to find the opposite point of the polygon where the line string of the 2 known points intersects the polygon.
Below is the SQL I'm using to find that point by taking the line string and extending it out past the polygon then getting the last point of intersection.
st_astext(st_endpoint(ST_Intersection(pg.poly_geog::geometry,
    (
    select
        st_makeline(a,
        st_translate(b,
        sin(az2) * len,
        cos(az2) * len))
    from
        (
        select
            a,b,
            ST_Azimuth(b,a) as az2,
            ST_Distance(a,b) + 0.00001 as len
        from
            (select
                st_makepoint(c1.longitude, c1.latitude)::geography::geometry as a,
                ST_Centroid(st_astext(pg.poly_geog))::geography::geometry as b )a )b))))

This seemed to work until I started running more tests and found that it only works if the linestring between the points is in a certain heading. If the linestring has a heading of right to left it works fine and I can get the first intersection point and the above calculation of the opposite intersection point. If the line string is heading in the oposite direction for another polygon then I'm getting the the point returned from the same side as the intitial intersection.

If I switch around st_endpoint to st_start point then it just changes around the direction that the function works for.
Ive also tried ST_GeometryN as well but that only works in a certain direction also.
How can I get this to function regardless of the heading on the linestring or how can I determine whether I need to use startpoint or endpoint


texts: AxTOCControl - Using ArcMap context menu

texts:
Using Ax*Control's in a custom dialog that runs within ArcMap.

When I right-click on a feature layer in the AxTOCControl, and go to display the context menu, is it possible to display and make usable:
1) The ArcMap feature layer context menu:
{BF643199-9062-11D2-AE71-080009EC732A}
esriCore.FeatureLayerContextMenu
or,
2) A custom menu, but with one or more of the existing commands e.g.:
{18DF94D9-0F8A-11D2-94B1-080009EEBECB}
esriCore.LayerContextMenuItems, subtype:9  (Layer_Properties).

The crux of the issue is that the existing menu items don't enable; presumably because I'm not giving the hook they expect, or specifying the selected layer in the way they expect.

I have seen the examples, but all show creating custom commands, with a custom solution to use the map control custom property to store the selected layer.  I've tried using that approach with the existing commands, but unsuprisingly it doesn't work.

If you can't use the existing ArcMap TOC context menu commands w/ the AxTOCControl, what is the recommended way of providing a similar context menu for the AxTOCControl - rebuilding from scratch in code, including the property sheets and pages ?


texts: How to implement a temperature map using leaflet/openlayers?

texts:
I am new on Leaflet/OpenLayers. I have been searching the internet for a way to create a temperature map like the one below, but I didn't succeed...the only thing I made was a classical heatmap, but not even closer to this example.

This is the scenario:

I have 5100 coordenates (longitude, latitude, metric)
This metric varies from 0 to 1

I want to implement a temperature map like the one below(from openlayers 3 cookbook):

I'm not looking for a 'ready to go' tool or a straightforward answer, I just want a 'direction'.

What I have implemented so far looks like this(completely different from the one above):


texts: Creating an aggregate function filtered by entry value

texts:
Can anyone help me parse an aggregate function for Qgis 2.18 which sums the values of one field, for all entries which match a  value of the entry in question.

So, for example if I have a simple layer with field "group", and "value", and entries are;

Group,Value

1,3

1,4

1,5

2,3

2,5

2,6

I want an expression which will yield "Group Total" such that,

Group,Value,Group Total

1,3,12

1,4,12

1,5,12

2,3,14

2,5,14

2,6,14

Next to all the entries in group 1 is the sum of all group 1, next to all the entries in group 2 is the sum of group 2

EDIT - I should add I've tried the most obvious sum("Value", "Group"), the predicted output in the expression calculator shows correct, but in the table the values are all just NULL, so I presume I'm doing something wrong.

EDIT2 - I should also have added that I'm looking for this function to automatically update somehow, in a virtual field, otherwise every time some value changes in my layer (and things do change!) all the totals will be wrong.


texts:
Given a distance in metres (x) and two coordinate pairs in lat/lon format (which form a line segment), I am trying to find a point x meters from point 1 on the line segment between the two points.

I know that in order to use Shapely for this operation, I need to transform the points. Does it make sense to convert the lat/lon pairs to UTM and then use interpolate? If so, will it work to use metres as the measurement for distance, or will I need to convert that to something else as well?

I'm open to suggestions on the best way to go about using Shapely for this.


texts: Create a mosaic from every available cloud free image over the study area and batch export from GEE

texts:
I have a study area which encompasses several images. I would like to create a mosaic for the entire study area and export it as a single image for every possible date for which imagery is available.
Here's what I'm using so far to create the mosaic. Most of this code has been taken from elsewhere and seems to fit the bill.
With the mosaics though, I want to calculate NDVI, threshold it and export just the thresholded NDVI band.
var batch = require('users/fitoprincipe/geetools:batch');

//add study area
var studyarea = ee.FeatureCollection('users/XXX/XXXXXX/studyarea');

Map.addLayer(studyarea, {color: 'red'}, 'Study Area');
Map.centerObject(studyarea);

// Filter the collection for the VV product from the descending track
/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

var Date_Start = ee.Date('2020-06-01');
var Date_End = ee.Date('2020-10-30');
var num_months_in_interval = 1;
var viz_params = {bands:"B4,B3,B2", min:0, max:0.3};
var images_to_display = 5;

// Create list of dates for time series
var n_months = Date_End.difference(Date_Start,'month').round();
var dates = ee.List.sequence(0, n_months, num_months_in_interval);
var make_datelist = function(n) {
  return Date_Start.advance(n,'month');
};
dates = dates.map(make_datelist);

function addSentinel2NDVI(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  //Thresholding
  //If NDVI less or equal to 0.3 => 0 else 1
  var thres = ndvi.gte(0.3).rename('thres');
  return image.addBands(mndwi).addBands(thres);
}

var CreateMosaic = function(d1) {
  var start = ee.Date(d1);
  var end = ee.Date(d1).advance(num_months_in_interval, 'month');
  var date_range = ee.DateRange(start, end);
  var name = start.format('YYYY-MM').cat(' to ').cat(end.format('YYYY-MM'));
  return ee.ImageCollection('COPERNICUS/S2_SR')
    .filterDate(date_range)
    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
    .filterBounds(studyarea)
    .map(maskS2clouds)
    .map(function(image){return image.clip(studyarea)})
    .map(addSentinel2mNDWI)
    .qualityMosaic('thres')
    .set({name: name});
};
var list_of_images = dates.map(CreateMosaic);

// Add a few monthly mosaics to the map.
for (var i=0; i<images_to_display; i++) {
  var month_mosaic = ee.Image(list_of_images.get(i));
  var label = ee.String(month_mosaic.get('name')).getInfo();
  Map.addLayer(month_mosaic, viz_params, label);

}

batch.Download.ImageCollection.toDrive(bandSubset, 'Mosaic',
                {scale: 10,
                 region: studyarea.getInfo()["coordinates"],
                 type: 'float'});

How do I specify that I want to export just the thres band for the mosaics as a batch process?
I know the syntax for it to apply to an image collection, unable to figure out in this instance.
//Export only the threshold band
//SYNTAX - imageCollection.select([old bands], [new band names])
var bandSubset = **VAR**.select(['thres'], ['thres']);


texts: Creating shortest route using line segment layer as network locations?

texts:
Using ArcGIS Desktop 10.2.1, I am trying to create a truck route that begins at one location on a county street grid and ends up in the same place, obeying a feature class of guard rail line segments throughout an entire county - we need the shortest possible.

I need to use these line segments as network locations. However, the only line-based NA objects that can be defined as network locations are line and polygon barriers. But I'm not trying to avoid each guardrail, just follow it.

I tried creating centroids and vertices from each segment to use as stops (on a single road as a test), but found out it sends me on this 1000 mile rat race to connect 3 points in a neighborhood road.

Is there a problem with how the tool is interacting with the LRS itself?

This seems like it would be a common problem, yet I haven't been able to find anything on it.


texts:
I have code that calculates the NDVI for agriculture plots.
now i'm trying to create raster that will be composed from 3 colores according to the NDVI values-
red= if the value is lowe then-> Mean- (Standard Deviation/2)
light green= if the value is between -> Mean- std/2 to Mean+std/2
Dark green= if the value is higher than mean+std/2

my code knows to calculate the statistics, my problem is the classification.
I tried to choose each time the pixels that are greater than those ''limits'' but for some reason, when I try to choose the light green values ( between Mean-Std/2 to Mean+std/2) it choose the pixels that are in the red category.

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
  .copyProperties(image, ['system:time_start']);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-06-01', '2019-06-30')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B1','B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

//test if clipping the image collection worked
Map.centerObject(geometry,9);
Map.addLayer(clippedCol.median(), rgbVis, 'RGB');

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

//analyze images from image collection collection
var listOfImages = withNDVI.toList(withNDVI.size());
var listOfNumbers =[0,1,2];

for (var i in listOfNumbers) {
  var image = ee.Image(listOfImages.get(listOfNumbers[i]));
  var meanDictionary = image.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry:geometry.geometry()
  });

 var STDDictionary = image.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry:geometry.geometry()

});

// print(i,'mean',meanDictionary,'standard deviation',STDDictionary);

var std2 = ee.Number(STDDictionary.get("NDVI")).divide(2);
var mean1 = ee.Number(meanDictionary.get("NDVI"));

var negBorder=mean1.subtract(std2);
var posBorder=mean1.add(std2);

var imageNDVI=image.select('NDVI');
var gtPOS=imageNDVI.gt(posBorder).selfMask().rename('PositiveBorder');
var ltNEG=imageNDVI.lt(negBorder).selfMask().rename('NegativeBorder');
var betMEAN=imageNDVI.gt(negBorder).lt(posBorder).selfMask().rename('MeanBorder');

var PositiveCOL = {
  palette: [
  '006622'
  ],
};

var NegativeCOL = {
  palette: [
  'ff0000'
  ],
};

var MeanCOL = {
  palette: [
  '00ff00'
  ],
};

Map.addLayer(gtPOS,PositiveCOL,i);
Map.addLayer(ltNEG,NegativeCOL,i);
Map.addLayer(betMEAN,MeanCOL,i);

 }

i'm also not sure if I can ''merge'' those 3 groups of pixels into 1 raster, so if you have any idea how to modify this code it will be great.


texts: Selecting PostGIS objects out of given distance

texts: I need to select by SQL querry objects that are further than 100m from specific type of road. For selecting object that are in given distance there is ST_DWithin but what to do in this case?


texts: Converting GK Northing / Easting into WGS 84 coordinates

texts:
I have a list of points in Excel where the header says GK Northing and Easting. I assume they are referring to the Gauss Kruger coordinate system.

Example:
GK Hochwert: 3497041
GK Rechtswert 5737303

How can I convert from this reference system to WGS 84? I have QGIS and Alteryx to perform this task.


texts: Calculating the maximum remoteness from of the boundary of a polygon (i.e. the radius of the largest circle inscribed in this polygon)

texts:
I have two related tasks. How is it possible to solve them (or at the at least the first one, which is the principle task) in ArcGIS or QGIS ?

Task 1.

How to find a radius of a largest circle inscribed in a polygonal object (shapefile) with a random shape of the boundaries? Under such circle means a circle that touches the polygonal object at two or more points. (Any circle with larger radius will intersect sides of the polygonal object).
Is it possible to solve this task in a vector?

If transform this task in a raster the question will be the following:

Task 2.

How to find for a group of cells of a raster (which represent the closed contour (a boundary) of a polygonal object from the first task which was put on a raster) a point (Latitude, Longitude) which lies at a maximum distance from all cells of this group and locates inside this closed contour? I need to determine the point’s coordinates and this “maximum distance”.

Are there any tools or scripts of this kind available in ArcGIS or QGIS?


texts: Distance from point to two neighboring points (both directions) in on polyline of pipe system

texts:
I need to find the distance of a valve to the two neighboring corners in a piping system.
I have the following shp files:
-Nodes for the valves
-Nodes for the corners
-line layer file of the piping system.
I need to compute the distance from a valve node, to the two neighboring corners. So the search should be in both directions. I looked at options to make a network of the piping system and find the distance from the valve to the two nearest neighbors, but the problem is the direction restriction.
If one searches for the nearest two neighbors the result could be two corners at the one end of the pipe that were the closed, instead of one corner at one side of the pipe and the other corner at the other side of the pipe.
I added a picture to illustrate the problem. I should get the distances from the valve to corners 1 and 2. But with the nearest neighbor I would get the distance to corner 1 and 3.

Does anyone know how to deal with this problem?


texts: Split Point Layer into even sized cluster

texts:
I want to cluster many points (around 1.000) into even sized groups (13 points each group). There is no attribute defining that groups, it should be grouped based on spatial distance.

I tried to work with "Clusterpy", "Concave Hull" and the long road via grids like in the following pictures. I made a four different grids (big to smaller ones) and the idea was to intersect the grid (starting with the big one) to intersect with the point layer. If there is the right number (13 points) or less but in the "surrounding area" not one point more as a single one - then that intersected point layer may be the first group. Then the progress of finding those spatial related groups would go on and the grid would get smaller and smaller. after the first grid i decided that this is not the ideal way to find a solution.

I also tried using "Shared Nearest Neigbhor Clustering" but with "13" or "8" as a number of Neighbors and only based on location every feature has the ID "1" in the resulting Output layer.

What would be the best way to split a point layer into even sized groups  (same feature amount) based on spatial proximity?

I use QGIS 2.18


texts:
Is it possible to incorporate an icon in the title toolbar of a dockwidget?

Semi-classification plugin has this. Is it possible?


texts: Suvey123 Saves Coordinates from Null Geopoint

texts:
I have a survey form with two geopoint questions. The first has bind::esri:fieldType set to null, and the second has the fieldType left blank. See a minimal example XLSForm below:

type
name
label
bind::esri:fieldType

geopoint
null_geopoint
Null Geopoint
null

geopoint
saved_geopoint
Saved Geopoint

When the survey is submitted, I would expect the coordinates saved as the feature's center, i.e. the point where the survey data shows up on a web map, to come from saved_geopoint. Instead, when I submit the survey, the feature is created at the coordinates of null_geopoint.
Switching the order of the questions fixes the problem, so it appears that the survey is always saving the coordinates of whichever geopoint question is first, rather than the expected behavior of saving the coordinates from the non-null geopoint question. In the full survey I am working on, through, it is not feasible to switch the order of the questions. Is there a way to get Survey123 to save coordinates of a geopoint question that is not first as its center?
(The survey is being published with Survey123 Connect version 3.11.123 and submitted with the field app version 3.11.164.)


texts: How to perform spatial join and obtain attribute in python/GDAL?

texts:
I am working on a script that will take an unprojected raster, and using a UTM grid (shapefile), will obtain the correct SRS (I have a shapefile where each zone is attributed with its numerical SRS) and perform a warp/projection using GDAL on that raster.  It is a quicker way than figuring out the right zone, and then manually doing the warp.

Procedurally, I am grabbing the raster corners and constructing WKT geometry (polygon); the objective is to do a spatial join with the correct feature in the UTM grid.  Here's what I have so far:

import os
import ogr
from osgeo import gdal

def getRasterBox(raster):
    src = gdal.Open(raster)
    ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform()
    lrx = ulx + (src.RasterXSize * xres)
    lry = uly + (src.RasterYSize * yres)
    return (uls, uly, lrx, lry)

def findUTMzone(raster_box, utmgrid):
    # do the stuff

utmzone_path = r'C:\temp'
utmzone_file = 'utm_grid.shp'
utmzone = os.path.join(utmzone_path, utmzone_file)

utm_ds = ogr.Open(utmzone)
utm_layer = utm_ds.GetLayer()

rasterdir = r'C:\temp'
raster_file = 'n38.tif'
raster = os.path.join(rasterdir, raster_file)

rasterdetails = getRasterBox(raster)
rasterWKT = 'POLYGON(({} {}, {} {}, {} {}, {} {}, {} {}))'.format(rasterdetails[0], rasterdetails[1], rasterdetails[0], rasterdetails[3], rasterdetails[2], rasterdetails[3], rasterdetails[2], rasterdetails[1], rasterdetails[0], rasterdetails[1])

Using the QuickWKT plugin for QGIS, I can see that I'm getting a good WKT ring for the raster bounds.  However, I'm unsure how to perform the spatial join in a way that grabs the feature attribute containing the correct SRS code.  How do I do this part (under def findUTMzone(raster_box, utmgrid):)?  I have seen a few arcpy examples, but am looking for something open-source.


texts: How to store selected attributes of vector layer using PyQGIS

texts:
I'm following this PyQgis workshop:
http://www.qgisworkshop.org/html/workshop/python_in_qgis_tutorial2.html#accessing-data-attributes

I have a point layer shapefile. I want to be able to retrieve and store x, y and name of each point( they already exist in attribute table and have indexes of 10, 11 and 1), hopefully in a dictionary format ((x,y) as key, name as value).

My problem starts here:

while provider.nextFeature(feat):
map = feat.attributeMap()

I can't figure out how to use this 'map' dictionary. iterating through it I either get

"PyQt4.QtCore.QVariant object at 0x23C9C3B0" or "built-in method toString of QVariant object at 0x10D38D50".

any help is appreciated!


texts: Building GDAL with libcurl support on Debian Linux

texts:
Does anyone know how to build GDAL with libcurl support on Debian Linux (specifically LinuxMint and/or Ubuntu)?

After following build instructions here and tip from @EvilGenius below, I'm specifying the path to libcurl with this command:

./configure --with-curl=-L/usr/lib/x86_64-linux-gnu

and I've also tried:

./configure --with-curl=/etc/ssl/certs/ca-certificates.crt

EDIT: I've built on a machine with PostGIS installed, for what it's worth, but I'm still not getting CartoDB as a data format in the list of supported formats.

Is there a way to check if GDAL was indeed built with libcurl support?


texts: QGIS 2.14.3 very slow when saving edits on Ubuntu 16.04

texts:
QGIS 2.14.3 "Essen" is really very slow when saving edits on a single layer (point features) on Ubuntu 16.04 64 bits.

There is only about 100 points in this layer. In editing mode, when saving, it takes about 15-20min (grayed QGIS windows during this time).
CPU (Athlon II X4 640) is always at 100% activity during this time.
Data are stored locally. There is no connection to any network.

What is the cause and how to fix it?


texts:
I am trying to make leaflet markers snap to a grid (e.g. 10px) when dragging.
The only method I can think of so far is to re-implement L.Draggable so that the offset calculated in _onMove is rounded to the nearest multiple of the grid size.

Is this the best way to do this? (Will it even work to override it then re-assign it to L.Draggable?)

I have seen plenty of questions about snapping features to another feature, which is relatively easy since you can just listen to dragend and do the snap at the end of the move. In this case I need to snap each 'drag'..


texts:
I have an ETL updating a layer attribute on a registered Enterprise geodatabase. The layer has been published on ArcGIS Portal which sends the layer with updated values to ArcGIS Online. This data is on a webmap that has Natural Breaks symbology applied to it. The web map participates on a web site. The layer is updated twice a day.
When the data values exceed the maximum value from the previous update the symbology can't render those now out of range values requiring someone to open the webmap and reapply the symbology each time the ETL runs.
Currently I'm using selenium to mimick a user logging on and doing the update manually but it is hacky and often errors out. Also requires an active user profile session on console machine.
I looked into Python API for ArcGIS and came across this:
https://community.esri.com/groups/arcgis-python-api/blog/2019/04/09/updating-layer-symbology-with-the-arcgis-api-for-python
but it looks like that only works if you already have the values for your breaks defined in an existing JSON file you use to update the current JSON file your web map is referencing.
Last, I tried to use arcpy to update symbology in an ArcPr project and export the layer file and get the new JSON values from that, but applying the code below did not update anything:
import arcpy

p = arcpy.mp.ArcGISProject(r'C:\my_project.aprx')
m = p.listMaps()[0]

for lyr in m.listLayers():
  if lyr.isFeatureLayer:
    sym = lyr.symbology
    if hasattr(sym, 'renderer'):
        sym.updateRenderer('GraduatedColorsRenderer')
        sym.renderer.classificationField = 'CaseRates'
        sym.renderer.classificationMethod = 'NaturalBreaks'
        sym.renderer.breakCount = 5

        lyr.symbology = sym

Seems like this should be way more simple than this, especially as this does not work.


texts: Seeking list of spatial functions available in PostGIS but unavailable in Spatialite?

texts:
It seems like most of the classic spatial function available in PostGIS (i.e. ST_Buffer) are usable with Spatialite, but some others don't seem to be supported in Spatialite (i.e.: no such functions: ST_ApproximateMedialAxis).

Is there a list available showing what is supported within QGIS DB-Manager for example?


texts:
I have downloaded the shapefile of Mexico's Admin Boundaries from the GADM site but I am not able to load it in QGIS.

It keeps giving me the following error,

  /Users/smaranh/Downloads/MEX_adm/MEX_adm2.prj is not a valid or
  recognized data source

But I simply downloaded the shapefile and loading it in the QGIS, then why is there a problem with the projection?

Can somebody guide me pls?


texts: Simultaneously exporting from multiple shapes by expression in QGIS?

texts:
I am working with OSM road shapes, there are 64 shapes for the whole American continent. I use "Select Features using an expression" tool to export selected classes from the "highway" field. (Expression is "highway" IN  ( 'motorway' , 'motorway_link' , 'primary' , 'primary_link' , 'trunk' , 'trunk_link' )).

Is there any better and quicker method to do this instead of manual routine?


texts: Custom Sensor type for New Ortho Mapping Workspace in ArcGIS Pro

texts: Is it possible to bring a sensor type (Satellite) which is not yet included in the available sensor list to create a new Ortho Mapping Workspace in ArcGIS Pro using Python ?? for example the sensor type for KhalifaSAT is not yet included in the list.


texts: I'm using R for the first time for a study on grasslands and I need to clip my findings onto a specific map. The trouble I'm having is the coordinates in the data don't match the coordinates on the map, they're mirror images. It's probably extremely simple since I can't find any tutorials or code online for it. Or I'm bad at wording this in search bars, either way, does anyone know how to do this transformation?


texts: How can I select only the points in a mixed geometry column in postgis postgres

texts:
Table "public.eating_drinking"
 Column |          Type           | Modifiers
--------+-------------------------+-----------
 osm_id | bigint                  |
 way    | geometry(Geometry,4326) |

I'd like to do something like select * from eating_drinking where is_point(way)='t';


texts:
Apparently if you stand at Greenwich then your GPS will not say you are at longitude zero. This is because the GPS reference meridian is about 100 metres to the east of the prime meridian.

My question is: Why is the GPS Reference Meridian 100m to the East of the Prime Meridian?


texts: Batch processing Vegetation Index with QGIS Raster Calculator

texts:
I'm applying the Green Leaf Index (GLI) to a multi-temporal set of RGB raster images to generate false colour overlays of the sections of the hay field. I can create GLI layers using the raster calculator and the RGB channels manually in the raster calculator. I'd like to speed up the process.
I'm looking for some pointers in how to set up a batch process in the raster calculator (or alternative process) so I can just enter the list of raster files and have the GLI generated for each raster on the input.
I'm just getting started on Python, so don't make any assumptions on the knowledge side.


texts: Keeping only ground points (class 2) with lasground in QGIS 2.14?

texts:
I want to keep only the LiDAR ground points in a las file, using the lasground tool from LAStools within QGIS. In QGIS Chugiak (2.4) there used to be an option to "keep class 2", but I don't know how to do it in the more recent version 2.14.

I tried to write in additional command line parameter(s) "keep class 2", but it breaks.

If anybody knows how to do it in lasground toolbox (bin) it would be also fine.


texts: Shapefile - problem with number of points in LinearRing

texts:
I am new to shapefiles and I've got a bit lost. I am trying to plot pooling sections in my region (Madrid) with different colors.

But I get

  Error: IllegalArgumentException: Invalid number of points in LinearRing
  found 3 - must be 0 or >= 4

when using the fortify function on my data (please, find the full code below). I found in this form that it may have to do with the rgeos library, but my rgeos version is

rgeos version: 0.3-8, (SVN revision 460)
GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
Polygon checking: TRUE

which is new, since I have installed it right now.

Someboty explained me that the problem is that my data has a point or a line section, not a full polygon, which requires 4 points. I went to QGIS, open the file and checked the geometry validity. I found the following errors

3969,segment 0 of ring 0 of polygon 1 intersects segment 0 of ring 8 of
  polygon 2 at -3.53569120482, 40.4258644642
3969,segment 0 of ring 221 of polygon 1 intersects segment 0 of ring 8 of
   polygon 2 at -3.53569120482, 40.4258644642
3969,Geometry has 2 errors.
4166,ring 1 with less than four points
4166,Geometry has 1 errors.

so that I elliminated polygon 1 and checked again the geometry validity of the new data set, but I get the same errors.

Can anybody provide me with a hint to fix this problem?

library(sp)
library(maptools)
library(rgdal)
library(ggmap)
library(ggplot2)
library(rgeos)
library(plyr)

seccionesOk <- readShapePoly('secciones_censales_buenas.shp')
seccionesOk@data$id <- as.vector(seccionesOk@data$DESBDT)

#Transform object of class "SpatialPolygonsDataFrame" in "data.frame"
seccionesOk.points = fortify(seccionesOk, region="id")

    ##### the previous line fails #######

seccionesOk.df = join(seccionesOk.points, seccionesOk@data, by="id")
colorsec=sample(1:4, size = length(seccionesOk@data$id), replace = TRUE )

cluster = data.frame(id=seccionesOk@data$id,
             seccionesOk@data$NAME_1,
             cluster=colorsec)

seccionesOk = merge(seccionesOk.df,cluster,by="id")

ggplot(seccionesOk) +
   geom_polygon(aes(x = long, y = lat, group = NAME_1, fill = cluster)) +
   scale_fill_manual(values =
     c("1" = "red","2"="yellow","3"="green","4"="blue")) +
      geom_path(aes(x=long,y=lat,group=NAME_1)) +
      coord_equal() +
      theme_bw() + xlab("Longitude") + ylab("Latitude")


texts: Seeking equivalent of ENVI's Change Detection Statistics in Google Earth Engine

texts:
I am looking to see LULC between two TIFF in ArcGIS Pro. In the class that I am teaching the students have access to pro 2.5, but this feature is only available in 2.6. I want to give them some exposure to Google Earth Engine.
I have 4 classes in each image that are equivalent between the two. I want to do a pairwise comparison to see percent change between each class pairing. This is easy in ENVI using the Change Detection Statistics tool, which the lab is based on.
If there is a way to achieve this in Google Earth Engine, what would the workflow be?


texts: How to install QGIS 1.8 plugin QGIS2KML?

texts:
I am trying to download QGIS2KML plugin.

I went to Plugins → Fetch Python Plugins in the Repositories tab and pasted the url.

It claimed that it was successfully downloaded, however, QGIS2KML was not installed in my  ~/qgis/python/plugins folder.

As a result from my failed attempt, I manually downloaded QGIS2KML and extracted it to my ~/qgis/python/plugins folder. Now QGIS2KML is shown as installed:

But I am still unable to see QGIS2KML icon.


texts: replaceDataSource() works in ArcMap but not as stand-alone script

texts:
I have one script (and one alternative) that I want to use to replace a broken data source. Neither is working outside of the ArcMap Python Window. Has it something to do with the MapDocument object?

Tested with a few different IDEs; including PyScripter, Pythonwin, PyCharm and IDLE.
In both instances the scripts run without error and print the correct print statements, but nothing is changed in the map document.
Only tested using Python 2.7.13 because when working with the arcpy module and ArcMap, Python 3.x is not supported. It's supported in ArcPro.
No difference with mxd.saveACopy() instead of mxd.save()

Script using replaceDataSource()

mxd = arcpy.mapping.MapDocument(r'\\gisfile\GISmaps\AtlasMaps\ATLAS_MAPS_18\Test
Folder\Beecher_Unit_200U_B.mxd')
for lyr in arcpy.mapping.ListLayers(mxd):
    if str(lyr.name) == "Test":
        lyr.replaceDataSource(r"\\gisfile\GISstaff\Jared\Hydro.gdb", "FILEGDB_WORKSPACE", "Hydro_WillCounty_Waterways")
        print(lyr.dataSource)
mxd.save()
del mxd

and using findAndReplaceWorkspacePaths()

mxd = arcpy.mapping.MapDocument(r'\\gisfile\GISmaps\AtlasMaps\ATLAS_MAPS_18\Test
Folder\Joliet_Elem_86_B.mxd')
mxd.findAndReplaceWorkspacePaths(r"\\gisfile\GISmaps\GISstaff\Jared\WillCoGIS_DataSHP.gdb", r"\\gisfile\GISstaff\Jared\Hydro.gdb")
mxd.save()
del mxd

On the other hand, these scripts work fine in the ArcMap Python window. For example,

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
if str(lyr.name) == "Test":
   lyr.replaceDataSource(r"\\gisfile\GISstaff\Jared\Hydro.gdb",
   "FILEGDB_WORKSPACE", "Hydro_WillCounty_Waterways")
   print lyr.dataSource
mxd.save()
del mxd

I've had a look at similar posts, but I couldn't find one helpful enough:

Same python code works in arcgis python window, but fails in pyScripter
Code Works in Python Window but not Script
Why does FeatureClassToFeatureClass work in ArcMap but not in Python script?


texts: Download zip from geoserver using javascript

texts: How can i download zip file, made by geoserver WFS request, implemented in javascript? i'm using OpenLayers.Request.POST to send request.


texts: Raster Resolution - Will not display past certain scale

texts: I just created a mosaic of a couple hundred digital elevation models in ArcGIS desktop. The mosaic only appears at the raster resolution. I am working with a very large area and need the scale of my map to be about 1:1,000,000. I am wondering what I need to do to have the raster appear in the area I need it to. I obtained the data from the USGS and it is already at a very fine resolution (1/3 arc second from the USGS, 1x1 cell in ArcGIS). All rasters were in the same coordinate system, same number of bands, etc.


texts: Using Python Script to get equal name as given and compare to multiple shp file in folder and copy as output for tool?

texts:
The following images will show my model. I used a piece of Python script and trying to get the output from it as a layer and use it for rest of the tool. However, I am getting error as displayed in one of the following image. could someone help me on solving this problem.

import arcpy
from arcpy import env
import os
arcpy.env.overwriteOutput = True
Match_FC=arcpy.GetParameterAsText(0)
Match_F=Match_FC[-13:-4]
#env.workspace =arcpy.GetParameterAsText(0)
output=arcpy.GetParameterAsText(1)
fcList = arcpy.ListFeatureClasses()
# Copy shapefiles to a folder
for fc in fcList:
    for fc in fcList:
        print fc[-13:-4]
        if fc[-13:-4]==Match_F:
            arcpy.MakeFeatureLayer_management('in_memory\fc', 'in_memory\output','')
            break

  ERROR 00229: Cannot open
  C:\User\Documents\ArcGis\Default.gdb\Rel_AA_001_Compari
  Failed to execute (Make Feature Layer)


texts: Programmatically loaded raster turns black after adjusting the georeferencing

texts:
I have an ArcGIS Add-In that scans our file directory for aerial imagery rasters intersecting a given point geometry and loads the qualified images to an mxd. As soon as you hit "Update Georeferencing" from the toolbar, the image turns black. If the file resource is loaded manually it renders correctly along with the updated georeferencing.

Before

After

Here are the two methods that deal with its symbology and loading:

private static void PullImages(IPoint centerpoint, ref List<IRasterLayer> rasterList, string directoryPath, ref int rasterCount, ref int imagesFound)
{
    IWorkspace ws; IDataset dataset; IRasterLayer rasterLayer;
    ws = RWSF.OpenFromFile(directoryPath, ArcMap.Application.hWnd);
    IEnumDataset enumDataset = ws.get_Datasets(esriDatasetType.esriDTRasterDataset);
    dataset = enumDataset.Next();

    while (dataset != null)
    {
        rasterCount++;
        OnUpdatePercentageLabel(rasterCount);
        IEnvelope rasterEnvelope = (dataset as IGeoDataset).Extent;
        rasterEnvelope.Project(NAD83);
        if (centerpoint.X > rasterEnvelope.XMin && centerpoint.X < rasterEnvelope.XMax &&
            centerpoint.Y > rasterEnvelope.YMin && centerpoint.Y < rasterEnvelope.YMax)
        {
            rasterLayer = new RasterLayerClass();
            rasterLayer.CreateFromDataset(dataset as IRasterDataset);
            RemoveHistogramStretchAndMinimize(ref rasterLayer);
            MXD.AddLayer(rasterLayer);
            Marshal.ReleaseComObject(rasterLayer);
            imagesFound++;
        }
        dataset = enumDataset.Next();
    }
}

internal static void RemoveHistogramStretchAndMinimize(ref IRasterLayer rasterLayer)
{
    IRasterRenderer renderer = new RasterStretchColorRampRendererClass();
    renderer.Raster = rasterLayer.Raster;
    renderer.Update();
    IRasterStretch2 stretch = (IRasterStretch2)renderer;
    stretch.StretchType = esriRasterStretchTypesEnum.esriRasterStretch_NONE;
    rasterLayer.Renderer = renderer;

    ILegendInfo legendInfo = (ILegendInfo)rasterLayer.Renderer;
    ILegendGroup legendGroup = (ILegendGroup)legendInfo.get_LegendGroup(0);
    legendGroup.Visible = false;
}


texts: why is the CDB_QuantileBins function not found in cartoDB?

texts:
What is wrong with the following SQL query on my data? I found pretty clear documentation (1,2) that the function is supposed to work.

SELECT CDB_QuantileBins(array_agg(numeric_column_of_mine), 10) from my_table


texts:
I am working with QGIS 3.2 (Bonn) and would like to add multiple pictures to a point (shapefile). In the end I would like to use the identify/info button and the picture should appear.
Anyone know how to do it without scripting? I am looking for an easy way to add pictures to the existing shapefile during the monitoring season, also multiple pictures


texts: Installing ArcPy (from ArcGIS Pro) with Anaconda?

texts:
I would like to use ArcPy from ArcGIS Pro with my Anaconda install.

ArcGIS Pro installs a 64bit python3 version of arcpy. How do I add this version of arcpy to my anaconda 64bit python3 install?

And if I need to find the arcpy version on the file system to do this, how do I find it? and how do I know that I've found the correct file?


texts: The grid labeling is not going on the border, it always marked on the boundary of vector dataset?


texts: Pickling a Classifier in Google Earth Engine

texts:
I'm busy trying to create a widely applicable reedbed vegetation classifier that could potentially classify any given area within Europe. I already have training data for an area around the Balaton lake in Hungary and the classification is very accurate >99%.

The problem is that the random forest classifier uses a lot of data and time while training. I know that within python there is a way to store a classifier using a "pickling" method, so if you would choose a different area you would not have to retrain your classifier every time. I was wondering if there is any way to do this in JavaScript within the Google Earth Engine?


texts: Divide all values in ASCII table (or raster file)

texts:
I have a pile of ASCII files representing 5min radar precipitation sums. The information is given in 1/10 mm which I need to convert into mm for further processing.
Is there a way of doing this, preferably as batch processing since I have more than 100 files? I tried converting the files to raster and use a con statement in the raster calculator (ArcGIS 10.5.1) which would not succeed...


texts: Automatically create donut polygons based on distance from boundary

texts:
I have a file with polygons representing risk of townships however the risk decreases inside each polygon so I'm wanting to create rings from each polygon for example from the boundary to 500m then another ring from 500m from the boundary to 750m from the boundary then another ring from 750m from the boundary to 1000m from the boundary and so on
I can do this manually by creating donut polygons but I was hoping for a more automated solution.


texts: Is aggregating points to areal units ever useful - Spatial statistics

texts:
I have a large point dataset of a population (By population, I mean that it was not a sample e.g. all houses in a county with an attribute such as size, price, age, etc.). Each point represents each case, their accurate position in space. It is possible to aggregate them up to any geographical unit such as census block, census tract, zip code etc. I would like to conduct several spatial statistical analyses (autocorrelation etc.) on this data. I tend to think that I should keep the data in the current format to maintain maximum variability rather than aggregating it to a spatial unit (e.g. census tract) to avoid MAUP and loss of variability.

However, is there ever a case to be made that under certain conditions a set of population point data should be aggregated to an areal unit?


texts: How to modify MODIS subset function to download product for a specific region and time period?

texts:
I am trying to assess the land cover change that has occurred over East Asia over last 18 years only in winters (Oct-Jan) for which I want to use MOD13Q1 dataset of modis from 2000-2018. I am using the mt_subset function to download the data.

I have already downloaded a sample data from Netherlands but I am unable to understand how to fetch my required dataset through this function which should have all the products from Oct-Jan (2000-2018) and the layer should be 250m NDVI. This website provides information about the library for downloading product https://cran.r-project.org/web/packages/MODISTools/vignettes/modistools-vignette.html

Here is how i am currently downloading the data.

VI <- mt_subset(product = "MOD13Q1",
                site_id = "nl_gelderland_loobos",
                band = "250m_16_days_NDVI",
                start = "2000-01-01",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = TRUE)

The location should ideally be as shown in the image (East China) and the start and end dates should only be for winters so that it only downloads products for that specific area and dates. How should I modify the arguments of this function to accomplish this task?