texts: Place image without using points in QGIS

texts:
I am trying to find a GIS software that does this, so I am trying QGIS but other suggestions are welcome.

I have a basic shp boundary such as a rectangle. I want to be able to import an unreferenced image (such as a tif or jpg) and place it as a new layer without needing to reference points immediately. I want to be able to use handles or another tool to re-size or stretch it in comparison to my shp. I realize this file is not necessarily referenced in this state even though it is floating over a referenced layer. It does not necessarily have to be saved as a referenced file ever as long as I can see it as an overlay while I work.

Basically because the image I am importing is not actual imagery and does not currently match anything on the map, I want to be able to place it in a way similar to placing and stretching an image in Microsoft Word.

Here's a sample image I photoshoped of what this feature might look like:

Is there a way to do this in QGIS, or suggestions of other software?


texts:
I am new to LiDAR point clouds and I have a big issue. I have LiDAR data for 30 km and I have 5 - 6 control points.

My point cloud is inaccurate (have averagely 3-4 m error in xy and 50 m error on Z) and I want to use the control points so that I can minimize this error. I have read about LAStool and FUSION, I just get information about shifting all point cloud in x,y,z but nothing about control points application.

Is there any software that can help me to apply control points?


texts: What is difference between GeometryService areasAndLengths method and Measurement widget in ArcGIS JavaScript API?

texts: What's the difference between the GeometryService areasAndLengths method and the Measurement widget to make a caculation in the ArcGIS JS API? Accuracy?


texts: How to change output folders when iterating a geoprocessing model?

texts:
I have digital elevation models clipped to watershed boundaries and I have to run a set of tools to hydrologically condition each one, and I want to use the iterator tool to automate this process.

Each DEM starts in a unique folder with the HUC 10 code, and each time I run the model, I want all the outputs to go into the unique folder that the DEM started in. Then, when the model iterates through to the next DEM, I want all of those output files to go in that unique folder. I have the recursive option in the iterator parameters checked so it will iterate through each subfolder, but I am having trouble figuring out how to set up the output so that each output file goes into the same folder that the original input came from.

Is there some sort of coding that I need to put in the output name for each tool? Something along the lines of %____%DEMfel.tif? I'm familiar with %n% and %scratchworkspace% but those don't work for me because it puts all the output files into one single folder.

I tried referring to this post, but that did not work for me: Iterating in ModelBuilder to Generate Output Files in Source Folder?


texts:
As the title: what is the most appropriate projection to reduce distortion, if I'm just showing the UK on a map?

Is WGS84 fine (since it's a relatively small area), or should I use some other kind of projection?


texts: Autopopulating field based on values from different field

texts:
When we ID our utilities (water pipes, storm manholes, etc.) we use ID's derived from a grid system. Our city is divided into grids that have a letter and number attached. "J14", "J15", "J16" and so on. That grid value is also assigned a number when a new feature is created. "J14-1" for example.
I already have a unique ID field that has the existing values. I would like to add to this when a new feature is created based on the location it falls within the grid system. We currently are switching asset management systems and no longer have the built in ID functionality that automatically did this. I have looked in to ModelBuilder and attribute rules for ArcGIS Pro but I have not seen a script that can do this and I don't have much Python experience.


texts: Structure shared raster data across GeoDjango apps?

texts:
I currently have a project that consists of multiple webmaps that are written as different apps in a GeoDjango project. The structure of the project is as such:
- project
  - manage.py
  - project
    - settings.py, urls.py, views.py, etc.
  - app1
    - templates, models.py, load.py, urls.py, views.py, etc.
  - app2
    - templates, models.py, load.py, urls.py, views.py, etc.
  - appX
    - templates, models.py, etc.

Each of the apps are focused on particular geographic regions, and are uniquely accessible to members from those regions (e.g., cities within a state) and I don't want to collapse them all into one app. That said, they all use the same raster datasets. At present, new raster tables are built for each app in the app-specific models.py file, and the raster data is split across these tables in the database.
I think it is more efficient to have a project level table that houses all of the raster data -- but am not sure how to do so.
So, rather than having app1_raster, app2_raster, appx_raster tables, I would just have a raster table that all three apps access.
How do I create a project level database table? From what I've seen, documentation encourages one not to do so.


texts: Given a flight path given in latitudes and longitudes, how do I find another flight path that is parallel to the path given?

texts: Given a flight route that is given in latitudes and longitudes, how do I find another route that is "parallel" to the route given (find a perpendicular line to each waypoint given)? I think the best way to do this is to transform the geographic coordinates to Cartesian and find the parallel route and then converted it back to geographic coordinates but I never get it right. I tried the following approach also but I am not 100% sure if this is correct and I never get the last point of the "parallel" flight path:


texts: How to clean vector segment intersect errors in GRASS

texts:
I'm using GRASS to process a vector map.
When I export back to .SHP and do a 'Check Geometry Validity' in Quantum GIS it lists several errors of this type...

3,segments 199 and 203 of line 0 intersect at 129334.076193, 669460.784787

I've tried using v.clean to get rid of them with no effect. Can anyone give me some hints for the correct v.clean tool to use please?

Also is there a vector geometry checker in GRASS so I dont have to keep exporting back to shapefile to check for errors?

Here are the commands I run. The errors start appearing in 'pass2' The map is free of errors prior to import into GRASS (as far as the Quantum GIS checker is concerned at least).

v.in.ogr -o dsn=E:/geography/uk2/gb_pass0.shp output=pass0
v.clean input=pass0 output=pass1 --o tool=rmarea,rmdac,rmdupl,rmline thresh=100000
v.generalize input=pass1 output=pass2 --o method=chaiken threshold=100
v.generalize input=pass2 output=pass3 --o type=boundary method=reduction threshold=100 -r
v.generalize input=pass3 output=pass4 --o type=boundary method=douglas threshold=40 -r
v.clean input=pass4 output=pass5 --o tool=prune thresh=40
v.clean input=pass5 output=pass6 --o type=area tool=rmarea thresh=7000000
v.clean input=pass6 output=pass7 --o tool=rmdangle,break,rmdupl thresh=1000
v.clean input=pass7 output=pass8 --o tool=rmbridge,rmsa,rmdac,bpol,rmdupl
v.clean input=pass8 output=pass9 --o tool=break,rmdupl,rmline
v.centroids input=pass9 output=pass10 --o
v.out.ogr input=pass10@mapset1 type=area layer=1 format=ESRI_Shapefile dsn=E:/geography/uk2/gb_boundary olayer=gb_boundary -e


texts: Is it possible to add a timeline using cartodb?

texts: I am trying to add a timeline in one of my visualization and I don't really know how to or even if you can add one in CartoDB. Does anyone know if you can and how?


texts: Problem using 'masker' while plotting KDE - R

texts:
The Problem

I am trying to plot the kernel density estimate for conflict data in Nigeria. I am using GISTools, rgdal, and sp loaded into R.

I want to generate the KDE, plot it out, and then mask the plot so that all excess space is clipped from the plot.

I am able to create the kernel density without problem. It's generating the mask that seems to be giving me the problem. I get the following error:

Error in `proj4string<-`(`*tmp*`, value = CRS(proj4string(input.poly))) :
Geographical CRS given to non-conformant data: -997.36173977 1014.70658626
-995.75375669 1013.91618459
In addition: Warning message:
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :
spgeom1 and spgeom2 have different proj4 strings

Code

Here is what I have done:

I import the shapefile for the administrative districts 2

nks_gis2 <- readOGR(dsn = "Map/Data/nks_admin2.shp")
nks_gis1 <- readOGR(dsn = "Map/Data/nks_admin1.shp")
nks_gis0 <- readOGR(dsn = "Map/Data/nks_admin0.shp")

I then check their CRS

proj4string(nks_gis2)

[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(nks_gis1)

[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(nks_gis0)

[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

I import the incident data in as a spatial data frame

nks_longlat  <- cbind(nks$longitude, nks$latitude)
nks_points <- SpatialPointsDataFrame(nks_longlat, nks, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

I check the CRS for the points

proj4string(nks_points)

[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

I now generate the kernel density estimate

nks_dens <- kde.points(subset(nks_points, country %in% "Nigeria"),
                   lims = subset(nks_gis2, NAME_0 %in% "Nigeria"))

And once again, I check the CRS

proj4string(nks_dens)

[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

I plot the density using level.plot

level.plot(nks_dens)

Now, here is where the problem occurs. I try to create a 'masker' so that it clips the plot to fit my Nigeria shapefile (I'll use the admin2 file):

nga_masker <- poly.outer(nks_dens, subset(nks_gis2, NAME_0 %in% "Nigeria"), extend = 1000)

I receive this error:

  Error in proj4string<-(*tmp*, value =
  CRS(proj4string(input.poly))) :

  Geographical CRS given to non-conformant data: -997.36173977
  1014.70658626 -995.75375669 1013.91618459

In addition: Warning message:

  In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
  unaryUnion_if_byid_false,  :   spgeom1 and spgeom2 have different
  proj4 strings

I have checked both my latitude and longitude, but do not see any long/lat that is even close to -997 or 1014.


texts: What kind of AWS configuration is needed for running a low-to-medium bandwidth Web map app?

texts:
Does anyone have experience with running Web maps (tile server + client JS scripting) on Amazon Web Services (S3, EC2 etc...)? What kind of AWS configuration is needed for running a low-to-medium bandwidth Web map app, covering a small(-ish) area (city to small-country size)?

All of the tiles would be pre-rendered and uploaded to S3. Ideally I would need a tile serving app on the web server that could serve MBTiles (instead of uploading hundreds of thousands of tile bitmaps individually). So some kind of EC2 instance would be needed, but what kind?

Thanks for any hints.

UPDATE: just to elaborate on my question. What I'm basically looking for is some feedback about how viable AWS is for hosting your own web maps as an individual (which means it shouldn't cost too much, say up to $30/month). I've been hosting my web maps for some time through "ordinary" hosting providers, but these have their own limitations (upload bandwidth is one, speed is another). I'm also looking for any good alternatives to AWS and anything to watch out for when using cloud services for web maps.


texts: Do FeatureLayers reproject their features to adapt them to the map projection?

texts:
I'm working with the ArcGIS Flex API.
For example, ArcGISDynamicMapLayers do reproject their data in order to make it in the same projection as the map.
Do FeatureLayers have the same behaviour?


texts: What are hardware requirements for a cloud hosted PostGIS & GeoServer?

texts:
I have been playing with a local installation of PostGIS & GeoServer. I would like to move my work online now to share it with several collaborators.

The proof of the concept map will consist of simple WMS and WFS services mapped through OpenLayers and GeoExt. It will be based on a map of ~3000 administrative units. Initially it needs to serve several people only. From then we will be working on polishing final view and functionality of the application.

There have been several suggestions made under answers to this question. For example, Nicklas is running his site on 512MB Linode (I think...).

I was thinking about using Ubuntu Server on cloud VPS, most likely hosted on rackspace.

My question: what would be the minimal server size I would need to use to start playing with this data online? There are few of them available. Once we decide to open the service to public rackspace promises seamless way of scaling up the service, so I would like to know what bare minimum RAM I need to start with.

Update 1: Similar discussion on GeoServer forum.

Update 2: Useful tutorials on setting up stack on AWS micro instance using Ubuntu and CentOS.


texts: Finding YMin and YMax of polygons in shapefile using ArcPy?

texts:
I'm using ArcMap 10.

I want a python script that  tells me the Xmin, Xmax, Ymin, and Ymax of a polygon. (I guess this is finding the extent of a polygon?)

My current script gives me the correct the Xmin and Xmax, but it gives me the Ymin and Ymax of the data frame. I want the Ymin and Ymax of a polygon.

Here is my current code:

# Specify the shapefile
shapeFile = "C:/Big_Creek_MapUnits.shp"

# Use the describe function
desc = arcpy.Describe(shapeFile)

# Print stuff out
arcpy.AddMessage(desc.extent.Xmin)
arcpy.AddMessage(desc.extent.Xmax)
arcpy.AddMessage(desc.extent.Ymin)
arcpy.AddMessage(desc.extent.Ymax)

Anyone have any solutions or tips?


texts: Is there a workaround to get WFS 'FeatureTypes' in a list with Python? I get timeout issues etc. with OWSLib

texts:
I have successfully been able to get the Feature type list from below answer, but only for smaller WFS datasets.

Importing all layers available on server for WFS in PyQGIS?

I need a workaround to get feature type list from larger WFS datasets due to eg. timeouts and loading time.
My current code is:

    from owslib.wfs import WebFeatureService
    wfs = WebFeatureService(url='https://.........../geoserver/ows?', version='1.0.0')
    list(wfs.contents)

Any ideas?


texts:
I am constructing a panel that includes past crop yields, but I found that SDA/WSS/GDG only grant access to the current version(2019) of cross-sectional tabular data.

Is it possible to access previous versions (dates)?


texts:
I'm trying to create an indoor trilateration program working off of N BLE objects, where I have any number of beacons in a room and my iPad can track its location within them to a sub meter degree of accuracy (I hope).

I asked this question before to check over my original Trilateration program. The only comment I got took me to this paper

Right now I'm trying to implement their iterative trilateration math in my code. Their formula is a follows (I am really sorry, I don't know how to put up the proper math formula):

  Reference points = Xi, Yi

  distance = di

  Trivial initial estimate = Xe, Ye

  Error in estimated distance =  |fi | =  di − Math.SquareRoot ((xi − xe
  )^2 + (yi − ye )^2)

Next they apply a delta to the initial estimate of Xe and Ye as follows:

  Xe = Xe + 0.05 Delta x

  Ye = Ye + 0.05 Delta y

With this new Delta, they then modified the original fi equation to look like this:

  |fi | =  di − Math.SquareRoot ((xi − xe )^2 + (yi − ye )^2) / di

Now, reading all that I've written my code as follows:

public GameObject refPoint1, refPoint2;
public float d = 0.0f;
public float f1 = 0.0f;
public float xe, ye = 0.0f;
public float multiplier = 0.05f;

// Use this for initialization
void Start () {

}

// Update is called once per frame
void Update ()
{
    xe = xe + multiplier * Math.Abs(refPoint1.transform.position.x);
    ye = ye + multiplier* Math.Abs(refPoint1.transform.position.y);
    d = Vector3.Distance(refPoint1.transform.position, refPoint2.transform.position);
    f1 = d - (Mathf.Sqrt((refPoint1.transform.position.x - xe) * (refPoint1.transform.position.x - xe)) +
                         (refPoint1.transform.position.y - ye) * (refPoint1.transform.position.y - ye)) / d;

}

In my unity scene I have placed the above code onto each of my three spheres twice so that each sphere is able to reference the other two. When I run this, however, my xe, ye and f1 just continue to grow and grow.

I'm just wondering if I've done something wrong, maybe misinterpreted something I shouldn't have? Any insight on this would be gratefully received.


texts: How to find the distance from a point to a polygon, via another point, in R

texts:
I am trying to calculate the distances between a central point, and the edge of a polygon, along straight lines drawn between the central point and another set of points (where necessary extended beyond the outer point to the polygon).

My data look like this:

library(sp)
library(raster)
# Centre point
Centre <- data.frame(x= 38.35419, y= 4.483533)
Centre <- SpatialPointsDataFrame(coords=Centre, data=Centre)
projection(Centre) <- CRS('+proj=longlat')
crs(Centre) <- "+proj=longlat +datum=WGS84"

# Sites
Site <- c(seq(1,10,1))
x <- c(38.29013,37.90987,38.29999,38.66006,38.15978,38.25978,38.15997,38.24062,38.25005,38.00712)
y <- c(5.33006,4.96017,5.23995,4.74010,5.09011,5.09015,5.06016,4.99723,5.06002,4.94634)
Sites <- data.frame(Site, x, y)
xy <- Sites[2:3]
Sites2 <- SpatialPointsDataFrame(coords=xy, data=Sites)
projection(Sites2) <- CRS('+proj=longlat')
crs(Sites2) <- "+proj=longlat +datum=WGS84"

# Polygon
Sr1 = Polygon(cbind(c(37.95169, 37.95169, 38.65000, 38.65000),
                    c(5.17153, 3.85601, 3.85601, 5.17153)))
Srs1 = Polygons(list(Sr1), "s1")
Hull = SpatialPolygons(list(Srs1), proj4string=CRS("+proj=longlat +datum=WGS84"))

I have managed to calculate the distance from the central point to each outlying point using:

spDistsN1(Sites2, Centre, longlat = TRUE)

Which returns the distance in kilometres.

To visualise the data, I use the code:

library(rgeos)
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, xlim=c(37.5,39.5), ylim=c(4,5), axes=T)
plot(Hull, col='green', axes=T, add=T)
points(Sites2$x, Sites2$y, pch=20)
plot(Centre, add=T, col="yellow", pch="+", cex=1.5)

and have managed to draw lines from the central point to each site:

C_P_lines <- list()
for (i in 1:length(Sites$x)) {
  C_P_lines[i] <- Lines(list(
    Line(rbind(Centre@data, Sites[i,2:3])))
    , ID=Site.coords[i,1])
}

C_P_spLines <- SpatialLines(C_P_lines, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

plot(C_P_spLines, add=T)

What I need to do, effectively, is to a) extend these lines to the edge of the polygon (for sites within the polygon), and then b) extract the intersection points of the lines with the polygon.

a) I have so far been unsuccessful in finding a solution to extending the lines.

b) I need to achieve the results shown here: http://uk.mathworks.com/help/map/ref/polyxpoly.html but in R. I have tried using gIntersection in the rgeos package,

Edge_points <- gIntersection(C_P_spLines, Hull, byid=T)
plot(Edge_points, col="red", add=T)

but this just crops the lines to the extent of the polygon.

I'm really new to using R for GIS work.


texts: QGis missing osgeo python-gdal, Python 27 vs Python 36

texts:
I have trouble while installing QGis from scratch on a computer that already have Python (3.4) + GDAL bindings. Here, QGis works with Python27, that needs its own GDAL bindings : how can i deal with this ? When i install a parallel version of GDAL, with a new GDAL_DRIVER_PATH, new Python27 paths etc... it breaks my Python36 install. It does that for any version of QGis...

Do you have any advice for that?


texts: How to generate multiple maps where one highlighted attribute changes?

texts:
I have a layer that contains counties within the state of Washington.

I need to build one map (jpg) for each county.  For that specific map all counties will show as light grey, but the specific county will be green, and the file name will be the county name.


texts:
How is it possible to add an image to a point in Qgis

I was thinking about something with html and actions

The following example is made by the plugin Qgis2web.
But the image wouldn't show. The image is stored in my documents;


texts: Avoid Overlap option not working in QGIS

texts:
I am creating polygons in a FGDB layer inside QGIS. I have topology on with advance options and have checked avoid overlap option for that layer.
Still, I am getting overlap some times and some times not when I am creating new features. Is this a bug or I am doing something wrong?


texts: Add Help Menu entry in QGIS 3 from `startup.py`

texts:
I would like to add a menu entry in the Help menu pointing to some web ressource, say https://gis.stackexchange.com. The following code executed from the python console works perfect:

from qgis.utils import iface
import webbrowser

def open_gis_se():
    webbrowser.open('https://gis.stackexchange.com')

iface.helpMenu().addSeparator()

gis_se_action = QAction('Go to gis.stackexchange')
iface.helpMenu().addAction(gis_se_action)
gis_se_action.triggered.connect(open_gis_se)

Result when typed in the python console:

... but putting it into my startup.py has no effect (Help menu remains 'as it is').

In QGIS 2, the above code put in the startup.py adds the desired menu entry as expected.

Why?


texts: Binary Raster (numpy array) to shapefile

texts:
i did a classification with a min/max threshold for stormdamaged areas. For that i calculated two indices (NSDVI) and (SAVI). I intersected two numpy arrays from those results. I recieve a binary GeoTIFF file.
For my end result i want to have a shapefile (GeoJSON or ESRI).
This is my code so far:

import numpy, sys
from osgeo import gdal
ds = gdal.Open('H:/Threshold/nsdvi.tif')
ds1 = gdal.Open('H:/Threshold/savi.tif')
rows = ds.RasterYSize
cols = ds.RasterXSize
myarray = numpy.array(ds.GetRasterBand(1).ReadAsArray())
myarray1 = numpy.array(ds1.GetRasterBand(1).ReadAsArray())
selection1 = numpy.logical_and(myarray1 < 0.4032,myarray1 > 0.0803)
selection = numpy.logical_and(myarray < 0.0005 ,myarray >-0.0255)
new2 = numpy.logical_and(selection==1,selection1==1)

driver = ds.GetDriver()
outDs = driver.Create('H:/Threshold/savi_nsdvi.tif', rows, cols, 1,     gdal.GDT_Float32)
outBand = outDs.GetRasterBand(1)
outBand.SetNoDataValue(0)
outBand.WriteArray(new2)

So i want to write all entries in the variable new2 which equal 1 as shapefile.

Afterwards i would like to filter out just those polygons>250m².


texts: GDAL VectorTranslate creates shards/fragments along anti-meridian

texts:
I have challenge to take a GeoJSON file and convert the CRS to WGS84/EPSG:4326 and also split the polygons along the anti-meridian.  These are both in preparation for presenting the data via mapbox.com.

I came across the VectorTranslate method of gdal which seems to miraculously achieve both of those requirements in 3 lines of code, e.g.

gdal.UseExceptions()
srcDS = gdal.OpenEx(datadir + 'input.geojson')
ds = gdal.VectorTranslate(datadir + 'output.geojson', srcDS=srcDS, format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])

My problem is the results are not quite perfect.  While it has split most of the polygons correctly, some are fragmented/shards. Here's screen shots of the input and output (via mapshaper.org):

Input:

Output:

The distortion was expected from the crs conversion, and looks ok on mapbox's background layers.  But some of the perfectly-aligned hexagons (from the input) which were split along the anti-meridian have fragments / shards missing, e.g. the highlighted polygon at the bottom of the output.

Does anyone know how to avoid or fix this?


texts: mapboxgl applying map.setFeatureState to a group of vector-tiles filtered features

texts:
Doing some tests, trying to let user click on a COUNTY and all counties belonging to STATE x get highlighted.

You can see running here

var map = new mapboxgl.Map({
        container: 'map',
        style: {
            "version": 8,
            "name": "blank",
            "sources": {
                "empty": {
                    "type": "vector",
                    "url": ""
                }
            },
            "layers": [{
                "id": "background",
                "type": "background",
                "paint": {
                    "background-color": "#1d1f20"
                }
            }]
        },
        center: [-82, 39],
        zoom: 3,
        debug: 2
    });
    var hoveredStateId = null;

    map.on('load', function() {
        map.addSource("counties", {
            "type": "vector",
            "url": "mapbox://reyemtm.28b260nr"
        });

        // The feature-state dependent fill-opacity expression will render the hover effect
        //  when a feature's hover state is set to true.
        map.addLayer({
            "id": "fill",
            "type": "fill",
            "source": "counties",
            "source-layer": "us_counties-8xelrv",
            "layout": {},
            "paint": {
                "fill-color": ["case", ["boolean", ["feature-state", "hover"], false], "whitesmoke", "firebrick"],
                "fill-opacity": 1
            }
        })

        function getUniqueFeatures(array, comparatorProperty) {

            var existingFeatureKeys = {};
            var uniqueFeatures = array.filter(function(el) {
                if (existingFeatureKeys[el.properties[comparatorProperty]]) {
                    return false;
                } else {
                    existingFeatureKeys[el.properties[comparatorProperty]] = true;
                    return true;
                }
            });

            return uniqueFeatures;
        }

        map.on("click", "fill", function(e) {
            if (e.features.length > 0) {

                var state_code = e.features[0].properties.STATE;
                console.info(state_code)
                var relatedFeatures = map.querySourceFeatures('counties', {
                    sourceLayer: 'us_counties-8xelrv',
                    filter: ['==', 'STATE', state_code]
                });

                var uniqueFeatures_state = getUniqueFeatures(relatedFeatures, "STATE");
                var uniqueFeatures_id = getUniqueFeatures(relatedFeatures, "id");
                console.log(uniqueFeatures_state)

                relatedFeatures.forEach(function(f) {
                    console.log(f.id) //18
                    console.log(f.properties.STATE) //allways the same, this is good
                    var s = {
                        source: 'counties',
                        sourceLayer: 'us_counties-8xelrv',
                        id: f.id
                    };
                    map.setFeatureState(s, {
                        hover: true
                    });
                });
            }
        });
    });

Console gives the correct info (that is, ids are unique and state_code is the same for all filtered features when click on one county). Curiously Texas is one of the few working correctly, but, when zooming in/out, highlighting becomes a chaos.

I have also tried using uniqueFeatures_state as the object to iterate over, with no success.


texts:
Possible Duplicate:
  How to create polygons with straight lines and right angles in QGIS?

In QGIS 1.8 is there a way to digitize buildings etc. with perfect 90-degree angles, Or, edit a polygon afterwards so that it has square sides?  I have a buildings shapefile.   I haven't been successful using the CAD Tools plugin.  Thanks to everyone with advice.


texts: Basemap annotations don't export to map image

texts: I am using the Open Street Map baselayer in QGIS 3.10 and when I export the maps to image in print layout, the basemap annotations don't print to the maps. The quality of the basemap is also quite poor and misses some of the details. Is there a way to fix this or can anyone recommend another baselayer to use?


texts: How to properly convert HDF5 file with 3K+ GCPs?

texts:
I am trying to create a GeoTiff from NOAA's Global RTOFS files for sea surface temperature, which have over 3000 GCPs, and any attempts at conversion using GDAL have resulted in slightly distorted output with strange artifacts at the sides of the image, as seen below:

What I have tried so far:
Using gdalwarp -t_srs EPSG:4326 HDF5:"rtofs_glo_2ds_f008_prog.nc"://sst out.tif, with various options added, such as different error thresholds and resampling methods, as well as forcing -tps (which resulted in even more distortion). All variations have pretty much resulted in the same output.
Using gdal_translate -a_srs EPSG:4326 HDF5:"rtofs_glo_2ds_f008_prog.nc"://sst out.tif, again with no visible difference in the output.
The output from gdalinfo -nogcp HDF5:"rtofs_glo_2ds_f008_prog.nc"://sst (the output was too long to post if I was to include GCPs) is the following:
Files: rtofs_glo_2ds_f008_prog.nc
Size is 4500, 3298
Metadata:
  Conventions=CF-1.0
  experiment=92.8
  history=archv2ncdf2d
  institution=National Centers for Environmental Prediction
  source=HYCOM archive file
  title=HYCOM ATLb2.00
  _NCProperties=version=1|netcdflibversion=4.5.0|hdf5libversion=1.10.1
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 3298.0)
Upper Right ( 4500.0,    0.0)
Lower Right ( 4500.0, 3298.0)
Center      ( 2250.0, 1649.0)
Band 1 Block=1125x825 Type=Float32, ColorInterp=Undefined
  NoData Value=1.2676506002282294e+30
  Metadata:
    sst_coordinates=Longitude Latitude Date
    sst_long_name= sea surf. temp.   [92.8H]
    sst_standard_name=sea_surface_temperature
    sst_units=degC
    sst_valid_range=-2.3368993 34.130985
    sst__FillValue=1.2676506e+30
    sst__Netcdf4Dimid=0

How are these kinds of files handled properly in regards to georeferencing?


texts: Converting an ENVI raster as point shapefile for editing in ArcGIS?

texts:
I have an ENVI raster with multiple classes inside. I tried to export this to a vector so that I can convert it to a shapefile in ArcGis (I want to manually edit/delete points in the raster), but the resulting .dbf file only has the following columns:

  OID | Class Name | Class_Id | Parts | Length | Area

It doesn't have x & y coordinates. How can I fix this and/or what are the other ways that I can export this ENVI raster to an editable point shapefile in ArcGis? Is this even the proper approach?


texts: Adding points along line by clicking on it using ArcGIS Desktop?

texts:
ArcGIS for Desktop 10.6.1.

I have a polyline (old railroad track). I'd like to add the train stops as points along it by clicking on the places. I have georeferenced historical maps as the background and can see the train stops on it.

The help files all point to the Construct Points tool, but with that tool, I have to select distances or equal distance, or number of points. I don't see an option for clicking on the line.

If I just generate points (create a new point file and just click), then then points may not be on the line.


texts: Is there a tool that performs automatic recognition of buildings?

texts:
I'm wondering if there is a tool, which

takes a satellite map of a region,
detects typical buildings and
allows the user to review its result (check, whether the tool correctly identified the buildings).

By "typical buildings" I mean the following: In some regions there are buildings, which look very similar from above.

For example, single-family houses have roofs, which look nearly identical (see red ellipse in the picture below).

Apartment houses also have a distinct shape.

So, theoretically it should be possible to identify houses and apartment buildings automatically (with human checking the results).

Is there a tool which allows this?

Update 1 (20.11.2013 10:47): Are there any LIDAR images for Russia, which are freely available or can be purchased at low cost?


texts: MapBox - one world copy on screen on high resolution screens

texts:
In my MapBox application, using MapBox.js Map is being displayed on full page (width:100%,height:100%). Map starts using following code snippet

L.mapbox.accessToken = accessToken;
map = L.mapbox.map('basicMap', 'mapbox.streets',
    {
        minZoom: 2,
        maxZoom: 17,
        attributionControl: false,
        worldCopyJump: true,
        zoomControl: false
    }).setView([35, 0], 2);

On large screens, World copies repeat itself (means USA and Australia comes twice from Left and right side). I Want to repeat world copies (worldCopyJump: true) BUT one copy at a time according to screens Like bounding box of [65,-100],[-65,100]

Tried options :
1) fitWorld : http://leafletjs.com/reference-1.2.0.html#map-fitworld

map.fitWorld().zoomIn();

2) fitBounds : http://leafletjs.com/reference-1.2.0.html#map-fitbounds

map.fitBounds([[
     65.958984,
     -100.353521
    ], [
     -65.50585,
     100.615985
    ]]);

But couldn't match the exact requirement. Any one has any idea, working solution for large and small desktop screens.


texts: Querying two Separate Columns in ArcPy SQL {where clause}?

texts:
Using ArcPy I'm Trying to Set up a Where Clause in my "FeatureClassToFeatureClass_conversion"

This {Where Clause} needs to Query Two Columns in the data:

FEATURE = 'ENVIRONMENT'

and

TEMPERATURE > 50

I'm having a hard time doing both queries in the same expression as you would using "AND"

  "\"FEATURE\" = 'ENVIRONMENT'" ?AND?  "\"TEMPERATURE\" > 50"


texts: How to get Not Installed Plugins in QGIS?

texts:
I am new to QGIS.

I am trying to export the Voronoi polygon Coordinates to CSV.

I have searched in net for that mmqgis plugin has to be installed but after installing the software in manage and plugin tools, not installed option is not showing.

How do I install that plugin and is there any other way to export?


texts: Is there an OSM Mapnik Gray Basemap in English

texts:
I have found a basemap, the OSM (Mapnik SW) map, that I would like to incorporate into a project in ArcMap. It has a great blend of gray scale methods to denote urban areas, terrain, and water features, while the gray scale leaves room for thematic layers.

The link to the map is here:

http://www.arcgis.com/home/item.html?id=19285b9249484a0585aa95e423b605da

My issue is that I need the map labels to be in English throughout the map. Unfortunately when I zoom into a non-English speaking region, like Taiwan, the basemap displays labels in Taiwanese.

Does anyone know of a similar basemap with English labels, preferably OSM Mapnik SW styled or a way to change these labels to English?


texts: Scaling NOAA CDR AVHRR NDVI imagery with image.reduceRegion()?

texts:
var reduce = img.reduceRegion(ee.Reducer.mean(), geom);

Using this code above on AVHRR NDVI imagery (with a resolution of 0.05 arc degrees), what scale should I set reduceRegion to? Scale accepts input of meters, but the resolution of this dataset is in arc degrees. Link to AVHRR NDVI metadata below:

https://developers.google.com/earth-engine/datasets/catalog/NOAA_CDR_AVHRR_NDVI_V4


texts: Getting the New Map When the FocusMapChanged Event Fires

texts:
Normally, in ArcMap, when an event handler to catch the FocusMapChanged event is set up, and a new map (data frame) is activated, when the event fires, a reference to the FocusMap returns the map that is being activated. BUT, when the currently active map is removed from the document, when the event fires, a reference to the FocusMap returns the map that is being removed! What?! This behaviour is not helpful to me.

Here's code that displays a message box with the name of the Focus Map whenever the FocusMapChanged event fires...

Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI

Public Class test_FocusMapChanged
    Inherits ESRI.ArcGIS.Desktop.AddIns.Button
    Private m_FocusMapChanged
    Private activeViewEventsPL As IActiveViewEvents_Event

    Private Sub ActiveViewEventsHandler(ByVal map As IMxDocument)
        activeViewEventsPL = TryCast(map.PageLayout, IActiveViewEvents_Event)

        m_FocusMapChanged = New IActiveViewEvents_FocusMapChangedEventHandler(AddressOf FocusMapChanged)
        AddHandler activeViewEventsPL.FocusMapChanged, m_FocusMapChanged
    End Sub
    Private Sub FocusMapChanged()
        MsgBox(My.ArcMap.Document.ActiveView.FocusMap.Name)
    End Sub

    Protected Overrides Sub OnClick()
        ActiveViewEventsHandler(My.ArcMap.Document)
    End Sub

Does anyone have any ideas of how to get around this and get a reference to the newly activated map rather than the removed map?


texts: Arcpy/Python how to change an attribute field based on another attribute variable

texts:
New to python here. I have looked through the 'similar questions' section to try and answer this myself and have found many attribute table python queries which are related to statistical operations. These typically call some other Arcgis tool. What I would like to do is a little different and I haven't found exactly how to do it yet.

I have two fields. The first has a list of name classifications for a cities general plan such as Lower Hillside, Residential Neighborhood, and so on. The second list houses a numerical 'code' which relates to conforming uses. Typically I have to select by attribute the specific general plan designation such as Lower Hillside and then use field calculator on the selected features to fill the second attribute table with the proper code.

What I would like to do is generate the codes based on the names. For instance, if the general plan column has "Lower Hillside" then the code column would equal "3" and so on. This would be repeated for a zoning designation as well. As pointed out in the comments this is a 1-to-1 generation. I could see something like the list of general plan names alphabetically assigned 0-x. I'm just not sure how to do that.


texts:
I have some map files consisting of 'polylines' (each line is just a list of vertices) representing tunnels, and I want to try and find the tunnel 'center line' (shown, roughly, in red below).

I've had some success in the past using Delaunay triangulation but I'd like to avoid that method as it does not (in general) allow for easy/frequent modification of my map data.

Any ideas on how I might be able to do this?

I'm working in fairly raw C++.


texts: Shortest Paths Between Adjacent Points on a Network

texts:
With the project I'm currently working on I need to calculate how much the total length of the network increases as I expand it from one node to the next. I have the road network as a line layer (32673 features) and nodes on the road network as a point layer (3461 features). A segment of these layers as viewed on QGIS is as follows:

As part of this project what I need to do is similar to v.net.allpairs in GRASS GIS but instead of a distance matrix connecting every combination of pairs on the network what I need is to connect pairs of adjacent (nearest) points along the network and find the length of each connection. I'm open to suggestions using any type of open source GIS software but prefer QGIS or GRASS GIS.


texts: QGIS displaying incorrect feature labels using perimeter

texts:
Using Qgis 3.16, I'm trying to correctly label buffers around an area with their radius in metres. I would like the labels to be located along each buffer's perimeter. All buffers are saved as a single layer, with radius as a parameter defined for each buffer. Therefore, the layer has labels set to 'Single Label/ Using Perimeter'.
The labels are displayed correctly when all features are visible on screen (see below).
However, when only a part of the map is visible (e.g. when creating a print composer layout), the labels get mixed/ superimposed on one another (see below).
This results in some features having double labels (the 5000m buffer, which is labelled with both '5000' and '7500' labels) and some having incorrect labels (all buffers >= 7500m, which are respectively labelled >='10000').
I have tried changing all labelling settings, including the ones under 'Repeating Labels', 'Priority' and 'Obstacles', but no joy.
Does anyone have an idea how this can be resolved?


texts: Shapely/ Python looping through a number of polygons

texts:
I am trying to loop through a set of polygons in shapely so using python) to determine whether they overlap or not. I do not want to do it manually as later i want to generalise the procedure and do this for a number of polygons imported from a shapefile (through Fiona).
These are the polygons I am now testing

polA = Polygon([(0,0), (3,0), (3,3), (0,3)])
polB = Polygon([(2,-1), (5,-1), (5,2), (2,2)])
polC = Polygon([(5,2), (8,2), (8,5), (5,5)])

I created a collection from these collection = [polA, polB, polC] and now I want a for loop to go through the list and check each polygon if it overlaps with the next. But nothing I try works, because I don't know how to loop and successfully test all the possible combinations. Don't know if I am understood. Does anyone have any idea how to do this using Shapely of Fiona?


texts: Layers not reappearing in composer after restarting QGIS

texts:
I would like to create several maps in QGIS, each showing various elements.  Layers and their styles are locked and in the same CRS.

Here is the rendering when I save the document before quitting:

And here is what I got upon reopening the document:

So how can I lock the OSM background? I'm using QGIS 2.18.7 and the layer uses the QuickMapServices plugin.


texts:
I am using the ESRI javascript API and I have a set of lat/long points I need to convert to a state plane coordinate system. My spatial reference is wkid 3418.

I can easily use the GeometryService to project the state plane to lat/long, but need to go in the opposite direction.

I was looking at the REST API for ArcGIS and saw a setting transformFoward. Is this the setting I need to go in that direction? (http://resources.arcgis.com/en/help/arcgis-rest-api/index.html#//02r3000000pv000000)

What are ways I can do this conversion?


texts:
I create my openlayers map like this:

tmpLayer = new ol.layer.Tile({
        source: new ol.source.XYZ({
            url: my_tile_url,
            crossOrigin: 'anonymous',
        })
});

The webmap I want to include only supports zoom levels 8-15. I want the map to disappear when the current zoom is not within that range.

I could find an example in the openlayers documentation where this is achieved with min/maxResolution. Is this also possible with specifying zoom levels?

To clarify, I do not want to constrain the zoom levels (like here), I want the map to disappear outside that range like in the above example.


texts: How can I get values from a array (which is evaluated by spatial query) as default-values in other fields

texts:
In QGIS 3.16.0 I have a postgis-polygon-layer with a field "array_polygon" which contains an array.
When digitizing a postgis-point-layer-point I would like to use the different array-values of the intersecting polygon in different fields as a default value to avoid lots of spatial queryies (one for every field).
The polygon layer has an "array_polygon"-field  with type "QVariantList" and type name "_float8" and length "-1" and precision "0".
The point layer has an "polygon_info"-fiel with type "QVariantList" and type name "_float8" and length "-1" and precision "0" too.
The "polygon_info"- field of the point-layer has a default value like
array_get(aggregate(
layer:= 'my_polygon_layer',
aggregate:='array_agg',
expression:=array_polygon,
filter:=intersects($geometry, geometry(@parent))),0)

The  preview shows array-values as expected:

But the "polygon-info"-field in my point-layer stays empty after digitizing a point over a polygon:

(I tried the expression without "array_get("...",0) too but without success).
My other fields - which should refer to "polygon_info"-array have a default value like
array_get( "polygon_info" ,0)

(I checked
array_get( current_value("polygon_info") ,0)

too)
but every approach has a NULL-result too.
What is wrong?

Has the expression or field-definition to be different to store the evaluated array (as seen in preview) in my point-attribute-table ?
Is the spatial-query not carried out before the default-values are set - how can I force QGIS to make the spatial query first and after that extracting the values from this field?
Or is my expression wrong - which alternative is correct ?


texts:
I am new to geo spacial analysis. I would like to learn kriging from a simple example. Say, we have four points with some z value.

library(sf)
library(gstat)
library(tidyverse)

dt<-tibble::tribble(
  ~id,  ~lon,  ~lat,    ~z,
  "A",   500,   500,  12,
  "B",  1000,   500,  13,
  "C",   500,  1000,  15,
  "D",  1000,  1000,  17
  )

Here I create an sf object and visualize it

DT_sf  <-  st_as_sf(dt, coords = c("lon", "lat"),
                 crs = 4326, agr = "constant")

ggplot() +
  geom_sf(data=DT_sf, aes(color=z), size=10)

How can I interpolate the z values, some average measurements onto the whole area?

So far I have explored I need to prepare a grid: like this (?) and a matrix of distances

grd <- st_sf(geom=st_make_grid(DT_sf), crs=4326)

dist<-spDists(as.matrix(dt[2:3]), longlat = TRUE)

coef = lm(log(z)~sqrt(dist), dt)$coef

I tried this and it is probably nonsense (i do not know what I am doing).

k = krige(log(z)~dist, as_Spatial(DT_sf), as_Spatial(grd), vgm(.6, "Sph", 900), beta = coef)

What package and command would result in extrapolation of z values onto the whole surface? What parameters would be needed?

The krige function now returns error:

Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  :
  NROW(locs) != NROW(X): this should not occur
In addition: Warning messages:
1: 'newdata' had 100 rows but variables found have 4 rows
2: 'newdata' had 100 rows but variables found have 4 rows


texts: Non-ASCII characters not displaying correctly after Python manipulation

texts:
I'm using QGIS 2.18 on a Mac (Sierra).

I'm doing some pretty extensive data manipulation using a virtual column, and using the function editor (so Python) to do so: one of my functions, for example, is to change street-name abbreviations to long-format (e.g.: All => Allée), and in French.

The result is encoded in UTF-8, and I can't seem to find a means to translate it into whatever encoding is needed to display it correctly on a map... The result of any string containing an accented character is 'null' in the attributes table.

Yet in another 'street name' layer, pulled directly from a PostGIS database, even accented labels (street names) display correctly.

Is this a common problem, and is there any way around this?

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from qgis.core import *
from qgis.gui import *
from qgis.utils import qgsfunction
import re

@qgsfunction(args='auto', group='Custom')
def streettype_format(input, feature, parent):
    rep = {
        "CITE":"Cité",
        "ALL":"Allée",
        "RUE":"Rue",
        "RPT":"Rond-point",
        "SENT":"Sentier",
        "VLA":"Villa",
        "TERR":"Terrace",
        "CHEM":"Chemin",
        "CAR":"Carrefour",
        "HAM":"Hameau",
        "BD":"Boulevard",
        "AV":"Avenue",
        "CRS":"Cours",
        "VOIE":"Voie",
        "CHAU":"Chaussée",
        "ARC":"Arcade",
        "GAL":"Galérie"
    }
    rep = dict((re.escape(k), v) for k, v in rep.iteritems())
    pattern = re.compile("|".join(rep.keys()))
    input = pattern.sub(lambda m: rep[re.escape(m.group(0))], input)
    return input

Also: from some d*cking around in the Python console:

>>>> import sys
>>>>sys.getdefaultencoding()
ascii
>>>>u"Héy"
u'H\xe9y'
>>>>"Héy"
'H\xc3\xa9y'


texts: Relation between Geoidal separation and antenna altitude in GGA sentence?

texts:
I have searched about definitions and relations about geoidal separation and antenna altitude reported in GGA sentence. The more I read the more I am confused.So in a simple manner can anyone explain me relation between altitude and geoidal separation ?

I am assumming that my gps source is using WG84 datum.


texts: Does Configurable Map Viewer (CMV) support custom symbology for operational layers?

texts:
I have looked through the documentation, browsed the widgets and searched here on GIS SE. I can't find any information regarding being able to override the default (published) symbology of a ArcGIS Server feature service.

Is this a feature that is available and I am overlooking it?

For now I am going to republish my service with new symbols, but like the idea of being able to alter it in the application to create different presentations of the same map.


texts: ESRI Web App - Publish publicly with login dialog

texts:
I have an ESRI web app hosted online I want to share outside my organization and the users on the other side don't have ESRI org handles. So I can't add them to a group and then grant them access.

I want to share this web page publicly, but to add security level of required user name and password via regular login dialog.

Is any other option within ESRI world before jumping into pure web development of implementing login dialog?
The plain vanilla way is to modify the generated JS code by adding login dialog. Eventually the app will be hosted online, but for the login part we need to set up local web server that will handle login requests, so we will have 2 servers accessed from single app. Is this something possible?

All in all seems like the thing I need is too popular use case.


texts:
I've generated a shapefile using .NET DotSpatial library, when I try to open the file in MapInfo I face below error, "Cannot open table. unable to open table road."

How do I determine what's the problem of generated dbf file?

When I open the dbf file using Open Office and overwrite the file after dummy change everything works fine.

Just uploaded the shapefile and dbf file here

Update

This is the helper class which used to generate shapefile.

public class ShapeFileWriter
{
    public static void Write(
        IEnumerable<IGeometry> geometries, DataTable dt, string path, bool overwrite = true)
    {
        var shapeFile = new FeatureSet
        {
            Projection = ProjectionInfo.FromProj4String("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +no_defs");
        };
        var columnNames = new List<string>();
        shapeFile.DataTable.TableName = "road";

        foreach (DataColumn column in dt.Columns)
        {
            if (column.ColumnName.Length > 11)
            {
                column.ColumnName = column.ColumnName.Substring(0, 9);
            }
            shapeFile.DataTable.Columns.Add(column.ColumnName, column.DataType);
            columnNames.Add(column.ColumnName);
        }
        int index = 0;
        foreach (IGeometry geometry in geometries)
        {
            IFeature addedFeature = shapeFile.AddFeature(geometry);

            addedFeature.DataRow.BeginEdit();
            foreach (string columnName in columnNames)
            {
                object o = dt.Rows[index][columnName];
                if (string.IsNullOrEmpty(o.ToString()))
                {
                    Type dataType = dt.Columns[columnName].DataType;
                    o = dataType.IsValueType ? Activator.CreateInstance(dataType) : null;
                }
                addedFeature.DataRow[columnName] = o;
            }
            addedFeature.DataRow.EndEdit();
            ++index;
        }

        shapeFile.SaveAs(path, true);
    }
}


texts: Similarity between two or more trajectories

texts:
I have the data of trucks (http://www.chorochronos.org/).

This data are gps coordinates of multiple trajectories of trucks in Athens.

I have to calculate the similarity between the trajetories, in order to delete those that are very similar!

Red And Green are similar, but blue, black and (red or green) are different trajectories.
I want to delete one of the similares, red or green.

Data are in points (geometry , lat and long ,  x and y)(coordinates gps), the image are examples of trajectories


texts:
I am doing a project where I am digitizing a map of Haiti from the 70s. Through R2V I digitized the lines and then I exported it into ArcGIS 10. I am just interested in the land use change, therefore, I have deleted the other vectors. I want all of the land use change areas to be complete polygons. I performed through data management in the toolbox a feature class to polygon operation. This was somewhat successful. However, it did not identify all polygons.

When I have complicated and nested polygons how do I get ArcGIS 10 to identify them as such?


texts: Can I use the globe plugin on a Windows OS and where can I find the plugin for QGIS?

texts: I am trying to add my layers (health facilities and pertinent infrastructure within a region) to a 3D image of the landscape and someone recommended I use the Globe plugin. I can not find the Globe plugin. Does anyone know where I can find it?/ Will the plugin work for a Windows OS?


texts: ArcGIS 10.1 crashes on use of python window

texts: When I open the python window in either ArcMap or ArcCatalog (ArcGIS 10.1) as soon as I press a key it freezes the entire program and it crashes.  It doesn't produce any warning, error, or pop up with info afterwards it just crashes.  I tried searching for similar problems through the forum but none of their solutions worked.  I installed Service Pack 1 and that didn't solve the problem.  I can use IDLE and/or PythonWin to import Arcpy but some custom python tools won’t work.  Some of those tools crash it in the same way.  Has anyone experienced anything like this or have any idea on how I can get python window to work for Arc 10.1? (I have also tried uninstalling and re-installing ArcGIS 10.1)


texts: GDAL install on Mac: pip doesn't see gdal.h

texts:
I'm trying to install GDAL on a Mac (macos 10.12.2). I downloaded the GDAL 2.1 Complete package from http://www.kyngchaos.com/software/frameworks (which installs GDAL 2.1.2), and successfully ran the installer. Now I'm trying to compile the Python bindings, using pip. In my bash_profile, I added lines for include paths (where the GDAL installer put them):

export CPLUS_INCLUDE_PATH=/Library/Frameworks/GDAL.framework/Versions/2.1/Headers
export C_INCLUDE_PATH==/Library/Frameworks/GDAL.framework/Versions/2.1/Headers
export LIBRARY_PATH=/Library/Frameworks/GDAL.framework/Versions/2.1/unix/lib

then I tried to compile the bindings with:

pip install gdal==2.1.2

but pip complained it could only find versions up to 2.1.0 ("No matching distribution found for gdal==2.1.2"). So I tried using 2.1.0 for the install:

pip install gdal==2.1.0

which got pip to attempt the install, but which eventually failed with:

fatal error: 'gdal.h' file not found
  #include "gdal.h"
       ^
 2 warnings and 1 error generated.
 error: command '/usr/bin/clang' failed with exit status 1

Yet, the file gdal.h is located in the Header directory I set for my include paths.

I tried the advice given in the following questions:

Python GDAL package missing header file when installing via pip

Install GDAL python package in virtualenv on mac

but it did not fix my problem.

What should I do next?


texts: PyQGIS, features geometry type changes to unknown after copying to another layer

texts:
I'm copying a subset of features from a shapefile based QgsVectorLayer (polyline type shapefile) to a memory based QgsVectorLayer('LineString', layername ,'memory'):

self.lyr_flightlines_selected.startEditing()

selected = []
for f in self.lyr_flightlines.getFeatures(QgsFeatureRequest().setFilterExpression(u'"id"={}'.format(selected_id))):
                feat = QgsFeature()
                feat.setGeometry(f.geometry())
                feat.setAttributes(f.attributes())
                selected.append(feat)

                print f.geometry().type(), feat.geometry().type()

self.lyr_flightlines_selected.dataProvider().addFeatures(selected)
self.lyr_flightlines_selected.commitChanges()

Selected layer features are presented on the map  without a problem, I can also access attributes of features I copied into the selection layer. But if I try to  access the geometry I get an  instance of QgsGeometry, but its type() returns unknown geometry. And if I print for example feat.geometry().asPolyline() I get an  empty list. I don't know why, because the features are visible on map?! hence they have to have valid geometry. What am I missing here?

(Using QGIS 2.8.2)


texts: Create common vertexes between shared polygon boundaries in PostGis?

texts:
Given two polygon tables L1 (blue) and L2 (red) that  share a common boundary but different amount of vertexes,   how can one make the shared boundary identical (they should share the same vertexes) ?

EDIT:

By identical I mean that 'If we say the Blue Layer is the reference layer, then the edge of the Red that is common with the blue should afterwards have 3 vertexes instead of two that is now'


texts: Wrap on Character is greyed out QGIS 2.12

texts:
In QGIS 2.12, my "Wrap on Character" option is greyed out. How do I make that active?
With it unavailable, I cannot get multi line labels.


texts:
I have a few scenarios where I need to dissolve vector data, and I'm wondering how to do this using PostGIS. What would the queries look like to perform these operations?

Below I'll describe the different cases through some simple examples.

The actual data I'm working with is considerably more massive and complex than these examples, in reality I'm working with data from OpenStreetMap (road networks, landuse polygons et.c.). PostGIS seems to be the foremost alternative for handling heavy data.

Segmented lines

Lines are split up into segments similarly to what is shown in the two images. The first line ("A-D") consists of three segments and has been split at B and C. I want to dissolve these three individual lines to get a line consisting of the vertices A,D. (vertices B,C are entirely redundant)

The second line isn't very different. It has been split at E and F into three smaller lines.

Dissolve polygons

The second issue concerns dissolving polygon data.  The geometries can either overlap or be within each other, touch (shared vertices) or nearly touch with some tiny distance separating the polygons.

For all of these cases I want to dissolve the polygons into larger geometries. I guess the last case where the polygons nearly touch is a bit special as the small space in between should be "filled out".

Final note

In the above cases the dissolving is solely based on spatial relationships. Non-geometric attributes aren't relevant.

As stated earlier the actual data I'm working with is pretty massive, and the base queries will likely have to be tweaked and optimized for handling very complex geometries.

If anyone would like to get the data shown in the examples it can be downloaded here as GeoJSON.


texts:
I am adding a picture marker symbol to graphic layer like this in Map click event:

const mapPoint:MapPoint = map.toMapFromStage(event.stageX, event.stageY);
mapPoint.spatialReference = map.spatialReference;
geom= mapPoint;

[Embed(source='assets/images/globe.png')]
var picEmbeddedClass:Class;

var pictureMarker:PictureMarkerSymbol = new PictureMarkerSymbol(picEmbeddedClass,100,100);

var myGraphicPic:Graphic = new Graphic(geom);
myGraphicPic.symbol = pictureMarker;
graphicsLayer.add(myGraphicPic);

which works fine. When I zoom in and zoom out it maintains its width=100 height=100 upon every change in zoom level to adjust its size.

How to make it retain its size as it was while adding it to graphics layer, irrespective of zoom level. As I zoom in it may seem bigger to me and small as I zoom out.
Thanks


texts:
I was attempting to follow a tutorial
to create a bubble plot of school district proficiency rates.

Data

shapefile of school districts (Census2010.shp)

District level data

Code

reading in data

dist <- read.csv('/Users/Path/to/District.csv', header = T)
state <- readShapePoly('/Users/Path/to/Census2010')

merge shapefile and district data

GA <- state
GA@data <-
    merge(GA@data, dist, by.x="NCESID", by.y="NCESID", all.x=T, sort=F )

creating bubble plot

plotvar <- GA@data$CCRPI
nclr <- 8
plotclr <- brewer.pal(nclr,"PuRd")
max.symbol.size=12
min.symbol.size=1
class <- classIntervals(plotvar, nclr, style="quantile") #Some school districts did not report data
colcode <- findColours(class, plotclr)
symbol.size <- ((plotvar-min(plotvar))/
    (max(plotvar)-min(plotvar))*(max.symbol.size - min.symbol.size) + min.symbol.size)
plot(GA)
GA.cntr <- coordinates(GA)
points(GA.cntr, pch=16, col = colcode, cex = symbol.size)
points(GA.cntr, cex = symbol.size)

Though the resulting map is just hollow school district outlines. What am I missing here?


texts:
Every time I run operations with temporary file output from the graphical modeler, I get the following error message:

  Oooops! The following output layers could not be open CLIP:
  C:\Users\Amine\AppData\Local\Temp\processing\371888fc2a8e444082652967444c577e\OUTPUTALG1.tif
  The above files could not be opened, which probably indicates that
  they were not correctly produced by the executed algorithm Checking
  the log information might help you see why those layers were not
  created as expected

This is a simple clipping operation that should not cause such drama in QGIS. I also get similar, if not worse errors, when outputting rasters using GRASS or SAGA commands. Can anyone please help? I've looked through and updated the processing options however  this doesn't change anything.


texts: Correlation between point feature (disease) and elevation

texts:
I have a problem! How can I calculate the correlation between point feature and Dem raster?

Hypothesis: disease prevalence in high altitude is more than low altitude.
I must do this in ArcGIS 10 and classify the altitude to 15 or 20 classes.

For example 0-200 meter above sea-level, 200-400 m and ... . I will join these altitude classes with disease prevalence shapefiles and then calculate the correlation between two entities. When I look at the layers, I see the correlation between altitude and disease.


texts:
I am trying to convert an .osm file to a .shp file.

I am using OGR2OGR.
There is an article at Openstreetmap which explains how this can be done. Here is an example command:

ogr2ogr.exe --config OSM_USE_CUSTOM_INDEXING NO -skipfailures -f "ESRI Shapefile" output.shp input.osm

And it works perfectly. As a result of upper command, in the "output.shp" folder I get four separate .shp files:
lines.shp, multilinestrings.shp, multipolygons.shp, points.shp

However what I would like to get is a single .shp file. So basically to convert an .osm file to a single .shp file.

Is it possible to do that directly through some change of the upper command?

If not then I tried to merge all these four .shp file into a single .shp file (merged.shp). I found an example topic here at gis.stackexchange. So I used the following 4 commands:

ogr2ogr.exe -f "ESRI Shapefile" merged.shp lines.shp
ogr2ogr.exe -skipfailures -f "ESRI Shapefile" -update -append merged.shp multilinestrings.shp
ogr2ogr.exe -skipfailures -f "ESRI Shapefile" -update -append merged.shp multipolygons.shp
ogr2ogr.exe -skipfailures -f "ESRI Shapefile" -update -append merged.shp points.shp

And the merged.shp file is created, but an error is raised at the upper last (fourth) command:

Attempt to write non-linestring (POINT) geometry to ARC type shapefile.

I tried to google but I couldn't solve this error.

Can somebody please help me?


texts: Export unique values from QGIS layer column

texts: I have a shapefile loaded into QGIS.When i apply filter i can quickly see all the unique values in that column. Is there a quick way to extract that for use in code? I know i could load the shapefile into a geo database and extract from there but for a single file that seems a lot of setup


texts: Variable-width (tapered) line styling in QGIS

texts:
I am looking for a way to style a line. I would like to change the line width between two vertexes, lets say from 2 to 4px. Is this possible to do in some way?

I was thinking about converting the lines to polygons but haven’t figured out any workflow that would work.

This is to make lines of different classes and widths melt together more beautifully.


texts: QGIS Modeler, making a date "Calender" choice parameter as input?

texts: Hello Currently I am trying to build a new model in the QGIS modeller. So far I have understand that I can predefine certain strings and with the refactor field give them as expressions in my new shapefiles. However I have one pecular question, for now there is a certain parameter where you have to type in a date "yyyy-mm-dddd".Of course you can just type the date in but I want to have a "Callander prompt" when you punt in date. How do you that?


texts: Determining optimum shared school bus route to a single destination from multiple origin points

texts:
Using ArcGIS Desktop and the Network Analyst extension, I want to determine the path of an optimized walking school bus route to school. a walking school bus is essentially a system by which adult volunteers chaperone elementary school children along a pre-set pedestrian path, to help keep them safe from traffic hazards and busy intersections. the idea would be that students would, starting from home, travel to their nearest walking bus "stop", then take the walking bus all the way to school.

I have a very detailed pedestrian network for my study area. scaled cost barriers include intersection controls, areas with high crime. the walking bus would bypass most of these barriers

Is there a workflow that could help determine the optimum location for the walking bus route itself? Right now my process is as follows:

Using closest facility analysis I determine the shortest path to the school for each house, then look at the density of routes along each segment of the network using a spatial join and the join count to find common "corridors". I would then digitize a route visually along this.

This is obviously fairly imprecise, and relies heavily on human judgement.

Is there a better way to go about this?


texts:
I am using ArcGIS viewer for flex v3.6 to create a web map. Can someone tell me if it is possible to print the overview map widget as part of your map?
P/s: I know someone asked this question 2 years ago but there wasn't an answer for that.


texts:
Every time I tried to export map with OSM map I get this error:

Traceback (most recent call last):
File "C:/Users/user/.qgis2/python/plugins\quick_map_services\py_tiled_layer\tilelayer.py", line 319, in draw
self.drawTilesOnTheFly(renderContext, self.tiles)
File "C:/Users/user/.qgis2/python/plugins\quick_map_services\py_tiled_layer\tilelayer.py", line 397, in drawTilesOnTheFly
tile_ds.SetProjection(str(transform.sourceCrs().toWkt()))
AttributeError: 'NoneType' object has no attribute 'SetProjection'

I export the map by the regular way- in the map composer. The  installsion process of QGIS 2.18 was just fine without any error message.

I don't understand what is my mistake?


texts: setting Arcgis world imagery as data source for a mosaic data set

texts: I have created a mosaic data set and want to set esri world imagery basemap as my data source but i am unable to connect to the server can anyone help


texts: How can I stop single pixel nodata values propagating through a tile pyramid when using gdal2tiles?

texts:
I am using gdal2tiles_parallel.py to generate a tile pyramid from a single geotiff with nodata value set. I have noticed that if I have a small hole in the original geotiff (say single pixel), then I get a single pixel hole in each higher zoom level I generate, so the hole effectively grows in geographic area as you zoom out.

I am using the average resampling method, and I believe the call which actually downsamples the code to generate each lower zoom level is gdal.RegenerateOverview; the gdalwarp documentation suggests that average resampling ignores nodata values, but I'm not sure that RegenerateOverview implements resampling in the same way and can't find any detailed documentation on how it works internally.

Does anyone know if RegenerateOverview ignores (or should ignore) nodata values, or if there is an alternative method I could substitute to handle nodata properly for resampling?


texts: How often are OSM data in Cesium updated with new OSM data?

texts: When users make edits to OpenStreetMap(OSM), or OSM data gets updated by any other means, how long does it generally take for those edits to show up in Cesium Ion or Cesium Sandcastle? Is this on the order of hours, days, or months? And what about other sources like Mapbox, ESRI data?


texts: Cannot define output with Concatenated path

texts:
I am working on a series of geoprocessing tasks and trying to define the output through the concatenation of folders, a geodatabase, and finally a raster dataset.  Each way I try it I receive an error stating: TypeError: unsupported operand type(s) for +: 'Result' and 'str'

Below is the code that I have so far with the Project Raster task giving me the error.

import arcpy
import os

#Set Environments
arcpy.env.scratchWorkspace = r"C:/GIS/temp/Subcatchment.gdb"
arcpy.env.overwriteOutput=True

#Create new folder and gdb as well clip and project DEM in WSPLN folder for new project

#Variables
Folder = "C:\GIS\TEMP\Subcatchments"
projectNbr = arcpy.GetParameterAsText(0)
gdbName = arcpy.GetParameterAsText(1)
watershedFC = arcpy.GetParameterAsText(2)
countyDEM = arcpy.GetParameterAsText(3)
outCS = arcpy.GetParameterAsText(4)
outTransformation = arcpy.GetParameterAsText(5)

#Create folder, gdb, and feature classes
projectFolder = arcpy.CreateFolder_management(Folder, str(projectNbr))
arcpy.CreateFolder_management(projectFolder, "MapDocs")
gdbFolder = arcpy.CreateFolder_management(projectFolder, "Geodatabases")
arcpy.CreateFolder_management(projectFolder, "Shapefiles")
gdb = arcpy.CreateFileGDB_management(gdbFolder, gdbName)

#Create buffer and clip county raster in memory
tempBuffer = arcpy.Buffer_analysis(watershedFC, "in_memory" + "\\wsBuffer", "2000 Feet", "FULL", "ROUND", "ALL")
tempRas = arcpy.Clip_management(countyDEM, "#", "in_memory" + "\\tempRas",
tempBuffer, "#", clipping_geometry="ClippingGeometry")

#Project raster to county coords and change units
DEMname = "DEM"
arcpy.ProjectRaster_management(tempRas, gdbFolder + "\\" + gdb + ".gdb" + "\\" + DEMname, outCS, "", "", outTransformation)


texts: How to highlight on mouseover without wiping out select control

texts:
I am trying to highlight features of a vector layer as the mouse moves over them. I can make this work.

I'm also trying to make the features selectable and have a popup box when they're clicked.
this works too.

What I can't make work is both of them at the same time. The hoverselect seems to override the popup select. I've based my code on http://dev.openlayers.org/releases/OpenLayers-2.11/examples/highlight-feature.html but with a few differences. Does anyone know why the hover select and select can't coexist?

select = new OpenLayers.Control.SelectFeature(vector_Layer, {clickout: true});
        vector_Layer.events.on({
            "featureselected": onFeatureSelect,
            "featureunselected": onFeatureUnselect
        });
        map.addControl(select);
        select.activate();

        select.handlers['feature'].stopDown = false;
        select.handlers['feature'].stopUp = false;

//REMOVED AS OVERRODE SELECT selectHover = new OpenLayers.Control.SelectFeature(vector_Layer, {hover: true});

function onPopupClose(event) {
    select.unselectAll();
}

function onFeatureSelect(event) {
        var feature = event.feature;
        var content = "
"+feature.attributes.SHEET.value + "
" + feature.attributes.description;
        if (content.search("


texts:
I'm trying to calculate the percentage of cell counts in new raster and the new raster is created from an old raster based on the Value > 14. In my code, 1st for loop is for original raster and keeping the value and count fields in a dictionary. Then, creating new raster using Value > 14. The 2nd for loop is for the new raster and keeping its value and counts of in a new dictionary and finally calculating the percentage.

However, I am getting an error on the 2nd for loop. The error says:

  'Raster' object is not iterable

Does anybody have a suggestion?

Code:

import arcpy, os
from arcpy import env
from arcpy.sa import *

#To overwrite output
arcpy.env.overwriteOutput = True

#Set environment settings
env.workspace = "C:/Subhasis/Test/Neshanic_Python"

outws="C:/Subhasis/Test/Neshanic_Python/extract"

#checkout ArcGIS spatial analyst extension license
arcpy.CheckOutExtension("Spatial")

# set local variable
inraster = ["01367620-r-r"]

for i in inraster:
    flds = ("VALUE", "COUNT")
    dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(i, flds)}
    sumcnt = sum(dct.values())
    newraster = ExtractByAttributes(str(i), "VALUE>=14")
    for j in newraster:
        flds1 = ("VALUE", "COUNT")
        dct1 = {row[0]:row[1] for row in arcpy.da.SearchCursor(j, flds)}
        sumcnt1 = sum(dct1.values())
        percentage=(sumcnt/sumcnt1)
        print percentage


texts: Modelbuilder variables to rename multiple shapefiles

texts:
I'm using ArcGIS 10.4 and created a model to rename a shapefile.

It works fine if the model is running from ArcMap but for only one time. Once tried to run it again, it gives error because the file already renamed from the previous run.

How can I make the shapefile name variable every time I run the model?


texts: SPOT 6 parameters to open data in matlab

texts:
I'm trying to open SPOT 6 data on Matlab using multibandread, but since I don't have some of the paremeters I'm only getting a black image.
I've been looking on the web for these parameters but I haven't found anything, not even in the User guide from Astrium. I have a tiff image.

The parameters I need are:

interleave (bip/ bil/ bsq)
Data type
byte order (le/be)
offset, xstart, ystart


texts: Merging large amount of polygon-features in script with GDAL/OGR Tools or Python?

texts:
We have a very large shapefile with over 30000 features and want to merge it with a script on a server. Currently we try it with QGIS, but the processing takes way too much time.

Is there a way to do it with the GDAL/OGR Tools or Python?


texts: Python code to determine the position on a map using Geo-Tagging

texts:
What is a proper python code, which enables to scan photographs (when captured from a series of snapshots) and match it with a geo-tagged image database.

Based on the extent of similarities existing between captured photo and geo-tagged images, it must help in finding my position on a map.

So far, I am aware of finding latitude and longitude of an image file using following python code

from osgeo import gdal

geotiff= gdal.Open('C:\Users\nikhil\Desktop\north_rgb.tif')

width,height=geotiff.RasterXSize,geotiff.RasterYSize

gt = geotiff.GetGeoTransform()

minx = gt[0]

miny = gt[3] + width*gt[4] + height*gt[5]

maxx = gt[0] + width*gt[1] + height*gt[2]

maxy = gt[3]

print(minx,miny,maxx,maxy)

I need some extension or modification of this code as my requirement stated above.


texts:
I have calculated Moran's I in R and reject the null hypothesis in both data sets, there is evidence of strong spatial autocorrelation.

However, how do I compare the two?

First set:

Moran I statistic
standard deviate = 9.7817, p-value < 2.2e-16

alternative hypothesis: greater

sample estimates:

Moran I statistic: 0.1736815445

Expectation: -0.0009970090

Variance: 0.0003188969

Second set:

Moran I statistic

standard deviate = 18.575, p-value < 2.2e-16

alternative hypothesis: greater

sample estimates:

Moran I statistic: 0.3307119647

Expectation: -0.0009970090

Variance: 0.0003189016

Both indicate a strong spatial autocorrelation, however, how can I precisely define between the two which one is stronger?

Is it simply by looking at which Moran's I value is greater?

Both have the same boundary (shapefile) and are in the same study area, but different time stamps.


texts: searching by property's value of geoJSON in OpenLayers

texts:
I have a map with 4 layers with the same attributes (not even the same values). an example of geoJSON is below:

{
"type": "FeatureCollection",

"features": [
{ "type": "Feature", "id": 0, "properties": { "CODIF": "PBA177432V1", "Position_X": 397632.884, "Position_Y": 419958.3456, "Point_Kilo": 177000, "Type": "2CmPt", "Voie": "V1", "dist_rail": 2.75, "Status": "PBA 3.50m" }, "geometry": { "type": "Point", "coordinates": [ -6.511578939607725, 34.374296518950914 ] } }
,
{ "type": "Feature", "id": 1, "properties": { "CODIF": "PBA177437V1", "Position_X": 397633.7158, "Position_Y": 419952.8751, "Point_Kilo": 177000, "Type": "STM.TT (PA)", "Voie": "V1", "dist_rail": 2.75, "Status": "PBA 3.50m" }, "geometry": { "type": "Point", "coordinates": [ -6.51156926129586, 34.37424727692602 ] } }
,
{ "type": "Feature", "id": 2, "properties": { "CODIF": "PBA177457V1", "Position_X": 397636.6028, "Position_Y": 419933.597, "Point_Kilo": 177000, "Type": "STM.TT (PA)", "Voie": "V1", "dist_rail": 2.75, "Status": "PBA 3.50m" }, "geometry": { "type": "Point", "coordinates": [ -6.511535636308449, 34.374073743241354 ] } }
]
}

I want to add a function that searches by value of property. ("CODIF" in this case). It means that when a user enters a value of property "CODIF" of a feature, he sees this feature.
I've added a function that searches by longitude and latitude values, but it has no relation with geoJSON properties (even if there are longitude and latitude properties in geoJSON file), as shown below:

$('#envoyer').click(function() {
    var lat = latitude.val();
    var long = longitude.val();

    if(long != '' && lat != '') {
        vectorSource.clear();
        vectorSource.addFeature(
            new ol.Feature({
                geometry: new ol.geom.Point([parseFloat(long), parseFloat(lat)]).transform('EPSG:4326', 'EPSG:3857')
            })
        );

        wkt.val('POINT(' + long + ' ' + lat + ')');
        map.getView().fit(vectorSource.getExtent(), map.getSize());
    }
    return false;
});

It responds to my need without being related to geoJSON properties values, and it's because of the nature of longitude and latitude properties, but for searching by "CODIF" property values, I must relate the function to geoJSON properties.

The HTML code of the function i should add is:

<div id="search">
  <input type="text" placeholder="Search PKs by CODIF">
  <button type="button">Search</button>
</div>

I started my function with something like:

$(document).on('click', '#search button', function() {
  var value = $('#search input').val();

What i need is to do a comparison of what a user entered as a value, with the values of "CODIF" property of geoJSON file, if it's found, the user will see the feature of this CODIF's value, and if not, he reads a message like: Not found.


texts: Extracting raster values and creating new raster in QGIS is not working as intended

texts:
I am trying to extract certain values from a raster image in order to create a new raster which should then only contain the specified values. I have a single band raster showing forest loss from 2000 to 2017. Pixel values range from 0 to 17, and correspond to the year in which the loss occurred. I am using the raster calculator with an expression such as:

"ForestLoss@1" = 11

It sort of works, but the output raster switches all the original pixel values to 1, and saves all of the other cells from the original raster as 0.

My hope was to have a new raster with only the cells that were originally 11 being retained.


texts: How to merge rasters with nodata cells in Arcgis

texts: I want to create a single attribute table base on multiple overlapping raster. Some raster do not overlap completely but have common areas. In that situation I need the value from other raster while keeping the non overlapping raster value as missing. I am using arcmap 10.


texts: Printing with Mapfish straight from Google Chrome without saving to file to computer?

texts:
I have Mapfish printing in ForestKit applicatioin. When I try to print a map, mapfish downloads the file to the downloads folder and then I can print the map.

It would be better for me if I could print straight from Google Chrome's printing preview.

Is this possible?


texts:
I am trying to give a polyline a new set of co ordinates from a server. When trying to use the .setLatLngs method to redraw the new line with updated co ordinates i get the error "Uncaught TypeError: a.setLatLngs is not a function". Am I using this method incorrectly? If so how do i use this method?

var a = new L.LatLng(0,0),
    b = new L.LatLng(0,0);

var polyline = new L.Polyline([a, b]).addTo(mymap);

//this is inside an ajax query
    a.setLatLngs([latA1,lonA1]);
    b.setLatLngs([latA6,lonA6]);

//latA1, lonA1 are new points for a point on the line
//latA6, lonA6 are new points for a point on the line


texts: 100 ArcSDE layers in one Map service - ArcGIS Server 10.3

texts: I am doing my very first WebGIS application with ArcGIS server using ArcGIS Javascript API 3.17. I have to bring 100 layers in my webgis application as REST service from Oracle DB. To do that whether can we publish 100 layers as one map service or each layer as one service i.e 100 map service or group the layers as 10 datasets and publish as 10 map service. please suggest that which is the efficient way to publish those 100 layers in ArcGIS server.


texts: How can I combine two vector layers with identical columns in virtual layer

texts:
I have two vector layers (shapefile) with the same columns and layout in each. One comes from the county and is frequently updated, and the other is manual additions I've made. I'm looking for a way to get a virtual combined layer for all the parcels from both layers.

It looks like QGIS virtual layer would work but I'm not sure how to use the query box. (very new to sqlite)

So far I have both layers imported into 'embedded layers' and I can get the data from one with:

select geometry as geometry, pin as pin from parcel_address

the other layer is named parcel_manual.

How do I get the data from both?


texts: domain for transverse mercator projection, using proj.4

texts:
I have to apologize for the naive question -- it's my first encounter with projections and this software package, and I intend to use it for a non-geographical application (hence I work on a unit sphere) -

when using the transverse mercator transformation, I find that it is not defined for latitude values between pi/2 and 3*pi/2, as well as -pi/2 and -3*pi/2 and thus discontinuous - why should that be, and how do I make this compatible with "my definition" of spherical coordinates, as (r,theta,phi) with -pi<=theta<=pi and -pi/2<=phi<=pi/2?
It's easy to see that for phi I can just substract pi/2 (thus set the equator to 0), but what about theta?

Thanks a lot,
Ben

import matplotlib.pyplot as plt
import numpy             as np

from pyproj import Proj

pi=np.pi

tmerc = Proj("+proj=tmerc +a=1")

def f(a,b): return tmerc(a,b,radians=True,errcheck=True)

nice   =[]
naughty=[]

for i in     np.arange(-2*pi,2*pi,.2):
    for j in np.arange(-2*pi,2*pi,.2):
        try:    a,b=f(i,j) ; nice.append((i,j))
        except: naughty.append((i,j))

plt.plot(np.array(nice   )[:,0],np.array(nice   )[:,1],color='green',marker='o',markeredgewidth=0.,linestyle='None')
plt.plot(np.array(naughty)[:,0],np.array(naughty)[:,1],color='red'  ,marker='o',markeredgewidth=0.,linestyle='None')


texts:
I am working on a query to find the closest point between polygon and point

Can I use this operator to calculate the distance:

SELECT a.point, b.poly FROM tableA a, tableB b
WHERE ST_Distance_Sphere( ST_Transform(a.point,4269), ST_Transform(b.poly,4269)) < 50
ORDER BY a.point<-> b.poly
LIMIT 1;

But it run and take a very long time without results or error message.


texts: QGIS: How to change symbology of a continuous raster to discrete with unequal intervals?

texts:
I'm new in QGIS and currently trying to convert my ArcGIS-layers into QGIS.

One of these layers is a slope gradient raster. I want to change the symbology so that a specific rang of gradients will be shown as one specific color in the map (pic on the right) instead of there being continuous colors (pic on the left).

I tried using the "singleband pseudocolor" rendering type, but the problem is that I only seem to be able to use equal intervals. Is there any way to use unequal intervals for my raster style?


texts: Apply symbols in a style sheet(stylex) through code

texts:
I'm in the process of converting an ArcMap add-in using ArcObjects to use ArcGIS Pro SDK. Not a fun process I must say.

We need to apply symbologies on the fly through code. In ArcObjects I can find the style sheet(.style file), loop through each symbol by name and return an ISymbol object which I can then apply to a ArcMap layer. Is there something similar in Pro SDK? For the life of me I can not find any example anywhere.

I was able to add the style sheet to a Pro project and apply symbols to a layer. With ArcObjects you can basically write code to do anything and everything but I'm not sure if this is the case with Pro SDK.


texts:
After exporting from QGIS2Web into the different formats available, I am going for Leaflet. I want to change the order of information displayed in the popup, to fix the position of the label above the marker and to change the boldness of the font (and if possible with a thin white line around the letters)
So far I managed to get the zoom I needed, to change the svg marker to the legend markers, to get a dashed line for some polygon.
Could someone point me into the direction or is there a description somewhere of what each file related to the index.html is doing ?


texts:
I want to use RasterCellValueCalculator in FME to perform an arithmetic operation on a pair of rasters (TIF).

FME requires that the two rasters must have the same number of rows and columns.

But my two raster have:

different extent
different size of pixel
different number of rows and colums

Which approach and transformers should I use to create a paired raster?


texts: How to aggregate multiple NC files into gridded data

texts:
I would like to convert daily total NO2 column data into an aggregate form similar to the monthly mean tropospheric NO2 columns on said website in order to later extract statistics using shapefile(s), but am currently unsure of how to best do this. Currently I'm trying to use harp to grid the data into L3 form, but my output doesn't look correct. Google Earth Engine provides a good example of functionality that I would like to achieve, except I want years earlier than 2018 (which is why I'm using the Temis NO2 dataset). For example, when running
harpmerge -a 'latitude>-90;latitude<90;longitude>-180;longitude<180' -ap 'bin_spatial(1801,-90,0.1,1801,-90,0.1)' QA4ECV* l3.nc
I get the strange image
Whereas a correct image would look like this:

And now I'm really stuck. Any tips and/or tricks, including libraries which aren't too hard to use?


texts: Referring to the parent layer with aggregate functions in QGIS

texts:
I'm trying to find how to refer to a field of a parent layer as a filter for an aggregate function. The instructions specify that the parent layer can be referred to with the variable @parent, but I've tried all manner of combinations and can't seem to get any to work.

The formula is;

aggregate(layer:='{my child layer}',aggregate:='sum',expression:="{a field from the child layer}",filter:="{a field from the parent layer}">1000)

Obviously, the field from the parent layer also appears in the child layer (otherwise the filter wouldn't work), but I want to build an expression which will iterate over several rows in the parent layer, each of which is a particular field in the child layer.

So far, I've tried replacing "{a field from the parent layer}" with;

@parent "{field}"
"{field}"@parent
"{field}@parent"
"{field}(@parent)"

...and probably a half-dozen other combinations of parentheses, quotes and orders.

The example simply uses 'geometry(@parent)', is it not possible with anything other than geometry?

Just in case the above doesn't make sense, or anyone can think of an alternative method, what I'm trying to do (in a simple example), would be if say one layer 'countries' had the fields "area", "urban population", "rural population", and "adult population". In each row is a different country.

In another layer 'summary', there are the fields "field", and "total areas". Each row is a field from the 'countries' layer.

I want the total areas column to add up the the total areas from 'countries' where those countries relevant populations are over 1000. So "total areas" for the row 'urban population' would show the total are of all countries where the urban population was over 1000, but "total areas" in the row 'rural population' would show the total area of all countries where the rural population was over 1000.

In the more complex real world problem I'm working on, there are 36 rows in the 'summary' table (36 fields to summarise from the child layer), and 9 measurements (in addition to area) to summarise. So writing a single query for each summary (360 of them) is just not an option.


texts:
I am trying to make a map of a terraformed Mars in QGIS and as a part of that, I need to convert the Mars MOLA DEM into contour lines at an interval of 100m. However, the file is too big for QGIS to process the entire planet so at best I have roughly 55% of the elevation data I need. I can switch to a 500m interval for the time being but what I need to do is somehow split my DEM in half - something I have not been able to figure out.

Does anyone know how I go about this?


texts: FME - FeatureReader parameters - Dataset: $(Output_GDB)

texts:
Novice FME user here.
I have an FME model that I need to run on a quarterly basis, which among other things writes an ESRI file geodatabase and calculates zonal statistics by way of a PythonCaller.
I've had to amend the model slightly to extract min and max values from the zonal statistics tables, where previously it only extracted the max values.

Currently models has a FeatureWriter (to fgdb) which passes to a PythonCaller (runs zonal statistics and saves results to tables in the fgdb) to a FeatureReader (extract results from tables to pass to AttributeManager).
Where I'm falling down is in the FeatureReader parameters, when I set the dataset parameter to $(Output_GDB) and then go to select the zonal statistics table outputs in "Feature Types to Read" I get a "Failed to retrieve feature types" error, essentially looking for an fgdb that doesn't exist.

I find that if I run the model, let it fail, then go back in and retrospectively point that FeatureReader to the fgdb created in the first run and assign the relevant values downstream, then delete the previous outputs and run it again, that it works.
But this is redundant and impractical.

Does the $(Output_GDB) value not look for outputs from upstream in the model? and How can I get around this issue?


texts: Creating Line (MultiLineString) layer out of CSV file

texts:
I need to replicate a *.qgs-file that has 4 different layers as shown in the picture below:

The layer that I am currently trying to replicate is the "links" layer. The information of the layer property is as below:

The attribute table of this layer looks like this:

The csv file, which is the source of the layer in question, looks exactly the same with the exception of the last two columns.
I was able to display the layer visualization using the "Points to Path", that comes from the "Processing Toolbox". However, I don't get the same attribute table as shown previously. In addition, it is not that straightforward as I had to manipulate the csv files beforehand.
I created a "To_From" column and created one csv file for "To_From" filled with "to" and one csv file for "To_From" filled with "from". After that, I add both files using "Add>Add Layer>Add Delimited Text Layer". Subsequently, I merge the layers and create a column "order" with 0 for "to" and 1 for "from". Only after this, I am able to use the "Points to Path".
I wonder if there is an alternative way for my current approach in which I am able to get the attribute table attached previously...
I am using QGIS 3.10.5-A Coruña.


texts: You set up your mxd in a folder.  However, when you copy and paste to a different folder, the Map Document Properties does not update the document path when you save a new mxd with a different title. It brings to the new folder from the previous folder ? I want to change that and how do I do that ?


texts: K-fold validation on Google Earth Engine

texts:
I have a satellite image "india_image" and data in form of a feature-collection "ft". I want to perform k-fold validation on this data in Earth Engine. How to do so?

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

var ft = ee.FeatureCollection('ft:1fWY4IyYiV-BA5HsAKi2V9LdoQgsbFtKK2BoQiHb0');

I am training a CART model using this data.

var training = india_image.sampleRegions(ft,['class'],30);
var trained = ee.Classifier.cart().train(training, 'class', bands);

I want to perform k-fold validation on this data. How can I achieve it?


texts: QGIS Select by expression no longer working as expected

texts:
I'm having a bit of trouble using "Select by expression" in QGIS.

I'm working on a shapefile of a polyline network, each line has a unique reference which is made up of letters, numbers and symbols such as 'CRCCU50234/50010'.

I often have to select lines(aka sections) based on a list of references so i've created a simple tool in excel to build a query to select by expression based on two inputs; a field and a value. The expression built looks like;

"Section_Re" = 'CRU50229/50010' OR
"Section_Re" = 'CRCC0229/50020' OR
"Section_Re" = 'CRCCU50229/50030' OR
"Section_Re" = 'CRCCU50230/50005' OR
"Section_Re" = 'CRCCU502/50011' OR
"Section_Re" = 'CRCCU501/50005' OR
"Section_Re" = 'CRCCU0232/010'

where "Section_Re" is the field name and the values in single quotes are the unique references.

I've been using this method of selecting sections based on a list for about 8 months now on multiple shapefiles without issue. A couple of days ago, this stopped working when used on any selection greater than about 10 records.

Strangely, the expression returns nothing when ran in full (for about 300 sections) but will return results if all but ~10 lines are removed from the expression. It's an OR expression so i would expect that even if some of the lines didn't return a result, those that did would be selected.

I'm very confused as to what's going on so I thought i'd turn to the stack exchange community for help. I hope that i've given enough information and explained my issue clearly, but if anything else is required, please let me know.


texts: Check is an asset (GPS route) has been on restricted road

texts: If I have GPS route for an asset and want to cross check the position of the asset against (pre-defined) "illegal" routes to have taken (series of GPS decimal), what might be a reasonable approach? My issue is similar to that of geofencing (triggers when an asset goes AWOL), but the restricted space is not a polygon representing forbidden territory, but multiple specific routes (in this case a route is the combination of a series of connected "illegal" GPS positions).


texts:
I have a .csv file with 3 fields: lat, long, and value.  I'd like to create a Thiessen Polygon raster from this file in Python without creating intermediate datasets.  Here's what I have:

# Create layer from points
arcpy.MakeXYEventLayer_management(points_file, 'lon', 'lat', "points_lyr", arcpy.SpatialReference(4326))

# Make polygons
arcpy.CreateThiessenPolygons_analysis("points_lyr", "tp", "ALL")

# Rasterize
PolygonToRaster_conversion("points_lyr", "grid_id", "weather_grid")

This fails at step 2, raising the following error:

arcgisscripting.ExecuteError: ERROR 000339: Input points_lyr does not have OIDs
Failed to execute (CreateThiessenPolygons).

Is it possible to create Thiessen polygons without first converting "points_lyr" into a shapefile?


texts: Searching for single Layer in Layers Panel of QGIS?

texts:
I'm working on a project that requires a lot of layers.

Because of this, even if i make different folders, it's sometimes hard to find the layers I want.

Is there a way to simply write a keyword somewhere to find (or highlight) the layer I want in the layers panel?

Is there a plugin for this?

For Example: I want to find the layer "Layer-Orange".

Where can I write "Orange" to find the specific layer?


texts: Changing Data Frame display units in ArcMap?

texts:
Last week, in ArcMap 10, I went to View -> Data Frame Properties... -> General tab to change the data frame map display units from Decimal Degrees to Degrees Decimal Minutes. (the Map Units were still in Decimal Degrees, in fact it is greyed out and can't be changed)

Today, I wanted to change the Display units back to Decimal Degrees, but after I click Apply and OK, the map status bar still shows coordinates in Degrees Decimal Minutes as I move my cursor around the map. I went back to Data Frame Properties and the Display units have reverted back to Degrees Decimal Minutes!

Why is this happening, and what should I do?


texts:
I use pyqgis in standalone app using python 2.7 in windows 10.

I have create some scripts using python idle and  work fine without some errors or messages.

Now I try to use pyqgis using pycharm editor with the some scripts where I have run it before in python editor.

in the pycharm that scripts work again,I get correct outputs but in the pycharm show me some errors messages and in the final show crashed error python but my output is correct.

like this :

QSqlDatabase: QSQLITE driver not loaded
ERROR: Opening of authentication db FAILED
QSqlDatabase: available drivers: QPSQL7 QPSQL QODBC3 QODBC
WARNING: Auth db query exec() FAILED
QSqlQuery::prepare: database not open
ERROR: Opening of authentication db FAILED

Process finished with exit code -1073741819 (0xC0000005) (and 255)

I don't get exit code 0 I get it -1073741819 (0xC0000005) and 255.

any idea why I take that error messages in pycharm and not in python editor ?