texts: How to align edges of overlapping polygons in the middle line?

texts:
The problem

Here are two polygons that overlap.

My ultimate goal is to reshape the grassland (green poly) and the reservoir (blue poly) so they share the edge somewhere in the middle of their intersection. It is acceptable if there would be a certain gap between the borders as long as they are not overlapping.

So, the area marked as W on the figure would be merged with the water body (blue poly to the left) and the area marked as G would be merged with the grassland (green poly to the right).

I have drawn the red line manually to illustrate the imaginary central line that splits the shared area into two equal (or roughly equal) pieces. The end result would look like this then (after merging):

Each polygon's shape was increased and they "are meeting in the middle".

I need to be able to do this analysis on any pair of overlapping polygons in an automated fashion. I am open to any existing solution. I am able to pull my data from an Esri geodatabase into other formats and do the processing somewhere else. If there are hints on the code libraries or tools for this kind of problem, I would be happy to learn more about them.

Research

There are some tools related to edge matching that exist in ArcGIS, mostly topology related. None of them supports splitting the overlapping polygons; they only allow choosing to what edge a polygon should be aligned to.

The Integrate tool would solve my problem in case the polygons are overlapping just over a minor area, but first the polygons can overlap over a significant area, and second the Integrate tool also modifies other vertices of polygons (not participating in the overlap) which I want to leave as is.

I got some crazy ideas about creating the longest central line over each polygon that would result when intersecting the water body and grassland, but I cannot really think through this using limited computational geometry knowledge I have. I guess the longest line that would split the polygon into two roughly equal parts would not always be the correct line I need to use.

I have started considering converting the polys to raster and reclassifying the overlapping pixels using an interpolation method (how far each pixel is from water and grass pixels), but can't figure out a workflow for doing this using ArcGIS Spatial Analyst tools.

Are there any good ideas on how to solve the problem?


texts: Developing shapefile (SHP) for location from scratch

texts:
I am working on a research problem for a location for which geospatial vector data does not exist. This means I have to create one from scratch so that I can create shapefile data for this location.
How can I go about this?
I would not mind texts and all kinds of tutorials. I am a physicist and I have not really worked with GIS data before, so consider me a beginner. I do not have any preliminary data yet but my location looks almost circular like this:


texts: How to keep empty entity during Geopandas Spatial Join

texts:
Is that possible to keep entities who haven't join entities in Geopandas Spatial Join ?

In

join = gpd.sjoin(buf, grouped_aquis, op='contains')

I would like to keep buf's entities that haven't grouped_aquis' entities joined.


texts: Calculate geometry breaks shapefile ArcMap?

texts:
Consistent error upon adding a field to a shapefile and calculating area. Area is calculated correctly and is displayed in the current ArcMap 10.3.1 instance.

Upon software reload, the shapefile cannot be added due to a mismatch between the features and data, and the file is effectively broken.

This problem happens locally and with shapefiles we have sent clients when they perform the operation on their machines. The shapefiles are created with mapexport from AutoCAD.


texts: Extract values from raster into points: data loss using extract function from raster package

texts:
I am trying to extract values of raster cells on points using extract from raster package in R (similar to 'Extract Values to Points' in ArcGIS-10.2). After doing so in order to check integrity, I computed the total value of all the raster cells (i.e. sum) ~66,000 and total extracted values on points ~49,000.  We have tried both techniques simple and bilinear interpolation techniques. Data loss gets only slightly reduced with bilinear technique, i.e. ~50,000 points are extracted.
There is a huge data loss, if anyone can provide a way forward or has experienced the same.


texts: Plot multiple Layers with mapview, Error with Gdal and crs

texts:
I have two rasters that I load with the raster-package. They look like this:

> ndvi
class      : RasterLayer
dimensions : 6547, 7496, 49076312  (nrow, ncol, ncell)
resolution : 0.0002899651, 0.0002899651  (x, y)
extent     : 19.96367, 22.13725, -29.91876, -28.02036  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

> slavi
class      : RasterLayer
dimensions : 6547, 7496, 49076312  (nrow, ncol, ncell)
resolution : 0.0002899651, 0.0002899651  (x, y)
extent     : 19.96367, 22.13725, -29.91876, -28.02036  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

When I call mapview(ndvi) + mapview(slavi) it gives me a superlong traceback.Besindes some warnings there is the error-message

Error in CPL_transform(x, crs, aoi, pipeline, reverse) :
  crs not found: is it missing?

and a little bit below

3: In CPL_crs_from_input(x) :
  GDAL Error 1: PROJ: proj_create_from_database: Cannot find proj.db

However, when I just plot one raster with e.g. mapview(ndvi) it plots it perfectly.

I am using gdal 3.0.4. I don't know if this might be an issue.


texts: How to trim strings in Field Calculator?

texts: I am  currently working on a County plat book and am in the process of prepping the "Owner_Nm" field that will serve as my annotation layer. I am trying to whittle down the owner names, based on parcel size, to either full name (easily done) to first initial + full last name or just initials. The format of the names is currently something like this "John H Smith" or "John H Smith & Lucy L Smith" for example. I know next to nothing about python in field calculator so I am hoping to get some help from this community. INRE: how to write my python code in field calculator to return, dependent on parcel size, either (using the above examples) "J H Smith" or "J H S" or "J H Smith & L L Smith" or "J H S & L L S". I have nothing to offer as to the code I've tried as it has mostly returned gibberish....Thanks in advance.


texts: How to re-project coordinates in a gml on-the-fly in Geoserver?

texts:
I am using Geoserver to get featureinfo, but sometimes I want to re-project the native SRS in the output gml.
The request is a (WMS) "GetFeatureInfo" with INFO_FORMAT: application/vnd.ogc.gml.

I am aware of the option to set the "declared SRS" when publishing a layer in Geoserver, see below:

But is there a way to re-project the coordinates on-the-fly?

EDIT: Below is two diffrent WMS GetFeatureInfo requests that shows the problem. The first request has Srs=EPSG:4326 and bbox in the same srs and the second request has in Srs=EPSG:900913 and the bbox in in the same srs.

http://demo.opengeo.org/geoserver/wms?REQUEST=GetFeatureInfo&EXCEPTIONS=application/vnd.ogc.se_xml
&BBOX=-139.848137,18.549615,51.852863,55.778385&X=351&Y=193
&INFO_FORMAT=application/vnd.ogc.gml&QUERY_LAYERS=topp:states
&FEATURE_COUNT=50&Layers=topp:states&Styles=&Srs=EPSG:4326
&WIDTH=780&HEIGHT=330&format=image/png
http://demo.opengeo.org/geoserver/wms?REQUEST=GetFeatureInfo&EXCEPTIONS=application/vnd.ogc.se_xml
&BBOX=-15567823.399227962,2101981.4669995294,5772234.305333375,7514424.280422591
&X=352&Y=220&INFO_FORMAT=application/vnd.ogc.gml&QUERY_LAYERS=topp:states
&FEATURE_COUNT=50&Layers=topp:states&Styles=&Srs=EPSG:900913
&WIDTH=780&HEIGHT=330&format=image/png

Both requests return the same gml with coordinates in EPSG:4326. I would expect the second request to have coordinates in EPSG:900913. As mentioned I can set the "declared SRS" in Geoserver, but what if I want to get the gml in another srs then the one declared. By on-the-fly I meant: Ask in one srs get the coordinates in the same srs without setting the declared srs. I hope that the question in more clear now.


texts: Coordinates of image corners for ImageExtent in Openlayers 5?

texts:
It is possible to get the coordinates of the four corners of the image, given the center of the image coordinate? or the bottom left and upper right corners?

Currently, I am center the map in the lat/long mentioned above and I want to display the image in that coordinates.

Example Code

import Map from 'ol/Map.js';
import View from 'ol/View.js';
import Projection from 'ol/proj/Projection.js';
import { getCenter } from 'ol/extent.js';
import ImageLayer from 'ol/layer/Image.js';
import Static from 'ol/source/ImageStatic.js';
import TileLayer from 'ol/layer/Tile'
import OSM from 'ol/source/OSM.js';

var osmSource = new OSM();
var extent = [0, 0, 4000, 3000];
var mapLayer = new TileLayer({
    source: osmSource,
    projection: new Projection(
      {
          code: 'EPSG:3857',
          units: 'm',
          axisOrientation: 'neu',
          global: false
      }),
      center : getCenter(extent)
  })

var map = new Map({
    layers: [
        mapLayer,

        new ImageLayer({
            source: new Static({
                attributions: '© dronemapper',
                url: './data/DJI_0014.JPG',
                // projection: projection,
                imageExtent: extent
            })
        })
    ],
    target: 'map-container',
    view: new View({
        // projection: projection,
        center: getCenter(extent),
        zoom: 2,
        maxZoom: 8
    }),
    projection: 'EPSG:3857'
});

As you can see the image is displayed in 0,0 (Guinea gulf), I'm trying to get the coordinates for bottom-left and right-up to display the image in the coordinates provided by the image.

PS: image info: 4000 x 3000 pixels
lat 38: 47: 40.98290000
lon 108: 1: 16.25819999

I have a working example with the image placed in the right coordinates, now I would need some directives to get the right zoom and rotation the image to match the road on the map.

Tip: the image if from here https://dronemapper.com ©, image DJI_0014.JPG

import Map from 'ol/Map.js';
import View from 'ol/View.js';
import TileLayer from 'ol/layer/Tile'
import OSM from 'ol/source/OSM.js';
import { Vector as VectorLayer } from 'ol/layer.js';
import { Icon, Style } from 'ol/style.js';
import VectorSource from 'ol/source/Vector.js';
import Feature from 'ol/Feature.js';
import Point from 'ol/geom/Point.js';
import { transform } from 'ol/proj';

function photoStyle(feature, scale) {
    var url = feature.get('url');
    return new Style({
        image: new Icon({
            scale: scale,
            src: url
        })
    });
}

function imgStyle(feature) {
    return [photoStyle(feature, 0.10)];
}

var layer = new TileLayer({
    source: new OSM()
});

var imgSource = new VectorSource();

var imgLayer = new VectorLayer({
   // extent: [0, 0, 4000, 3000],
    source: imgSource,
    style: imgStyle
    // ,
    // rotation: Math.PI / 6
});

//Photo coordinates   38.79471747, -108.02111111 var extent = [0, 0, 4000, 3000];

var feature = new Feature();
feature.set('url', './data/DJI_0014.JPG');
var coordinate = transform([-108.02111111, 38.79471747], 'EPSG:4326', 'EPSG:3857');
var geometry = new Point(coordinate);
feature.setGeometry(geometry);
imgSource.addFeature(feature)

var map = new Map({

    target: 'map-container',
    view: new View({

        center: coordinate,
        zoom: 19

    })

});

map.addLayer(layer);
map.addLayer(imgLayer);

You Can use [ALT] + [Shift] + Click to rotate the map.
I would like to display the image rotated and zoomed to the right position.

What I want

What I have


texts: How to export a map from QGIS 3.10.2 as a georeferenced map to be used in Avenza Maps?

texts: I am familiar with Esri ArcMap and ArcPro, but I am extremely green when it comes to QGIS.  The version I am using is QGIS 3.10.2 on macOS High Sierra 10.13.6.  I have created a map using google maps satellite imagery for forest inventory purposes. Now I wish to export it as a georeferenced map to be used on Avenza Maps.  When I go to export the map, the GeoPDF option is greyed out. So my question is why is that, am I potentially doing something wrong, or is there a better way to go about it?


texts: How to clip a points layer using a polygon/defined borders?

texts:
I have a global ocean temperature point data downloaded from NOAA.

I want to clip and isolate data just for Southeast Asia in the form of a rectangle.

I would also want to do this to other parameters such as salinity, phosphate, etc, so I would need to know the boundaries to be used for clipping the other parameters.


texts: Improving polyutil islocationonedge to return coordinates instead of true/false value

texts: I want to modify the PolyUtil.isLocationOnEdgeOrpath to return location coordinates instead of true/false value. I have tried different hack to use the returned value for another query whereby if response is true, use that coordinates that returns true to select all the coordinates from a MongoDB collection.


texts: Overlaying polygon on other polygons and creating new one using ArcGIS Desktop?

texts:
I would like to write a seeding plan. I created 3 zones and I want to insert a control strip (now I have 3 polygons, in the end I would like to have four). But I just can't merge the two together.

If I merge them, the program merges every polygon into one big one. I tried to run through the forums but I think I just can't write the question correctly.


texts: PNG from GeoServer size on Google Earth seems uncontrollable

texts:
I have set up a GeoServer instance that hosts a lot of data that we process at work. It seems to be the easiest way to distribute since I can provide a simple network link to the GeoServer so that others can view it in Google Earth. In this case, I have a bunch of points that I have created a style for. The style is a PNG for each type. Here is an excerpt from my SLD:

<PointSymbolizer>
    <Graphic>
        <ExternalGraphic>
          <OnlineResource xlink:type="simple" xlink:href="towerBlue.png" />
          <Format>image/png</Format>
        </ExternalGraphic>
        <Size>10</Size>
    </Graphic>
</PointSymbolizer>

The sizes of these seem to be all out of whack. When zoomed out, they are tiny, but when I zoom in literally 2 wheel clicks from where the first screenshot comes from, they cover the whole screen. What could I be doing wrong here? Also, SVG didn't work for me at all, which is why I'm using PNG.


texts:
I have two rasters, one with land use data (raster A) classified in discrete values that represent different classes (1,2,3,4,5...18). The other one (raster B) contains data that can represent precipitation, temperature, humidity, etc. In both rasters, each pixel has a different value.
Both rasters have the same coordinate system and resolution, although the extent varies. The first raster covers only a patch of the second raster.
I want to intercept both rasters, and where a certain class from raster A with a certain value from raster B intercept, it multiply raster B times a constant value that depends on A. I have these constant values in a CSV file.
For example, if raster A land-use category==7 intercepts raster B where the value of B is equal to 50, then multiply 50 times 0.15. Another example would be, where
raster A == 1
raster B <  50

then
raster C == raster B * 0.46

but, wherever
raster A == 1
raster B > 50 & raster B < 150

raster C == raster B * 0.30

Important to mention is that for the same land-use class, different default values exist that depend on the value range from raster B (as shown in the second part above).
The question is, how can implement such conditional operations in R? I have tried multiplying raster A so values are scaled up (x100) and I can freely manipulate the resulting merge, without a success.


texts: Restart autoincremental numbering by unique values

texts:
I'm using the graphic modeler and 'Add unique value index field' algorithm to assign IDs to unique records in a vector layer which works great, but then I want to add a 2nd field that will autoincrement for each unique value created.  I can use the 'Add autoincremental field' algorithm, but then I get 1 to n for my IDs, where I'd like it to re-start for each unique id.

Example: rather than 10-0, 10-1, 10-2, 11-3, 11-4  I'd like to see 10-0, 10-1, 10-2, 11-0, 11-1.

Bonus: If there's a way to start the autoincremental number at 10 instead of 0, that would be even better!


texts: GDAL Plugin with a different version of GDAL aready built into QGIS

texts:
My qgis windows version is 1.8.0 Lisboa and use the gdaltools in version1.9.2.

I want to use the plugin NetCDFBrowser, but this plugin need the gdal 1.10.

I install 1.10 stand-alone gdal found here http://www.gisinternals.com/sdk/
but i don't know how to make it works with qgis?


texts: I got 2 landsat 7 images from the same place but different date of adquisition, image A is well placed, but image B is shifted about 1500 meters east and 14700 meters north. I need to use both images to correct stripping issues, because of that both need to be aligned. How can i correct image B to be aligned with image A in grass 7?


texts: Which geostatistical tools to use to determine surface changes over 40 years?

texts:
I have a data from the Weather station with SWEmax (snow water equivalent) collected over 40 years. I need to interpolate the data over my river's watershed and then, somehow, determine the differences which might have happened over this time.

Which tools do I need to use in ArcGIS 10.1 to achieve this?
Do I need to make z-surfaces first and then take statistics from them?
Is the best way to show the differences is to use a diagram with years on the x-axis?

Thanks.


texts: Create a custom CRS based on a known PSM

texts:
We have some CAD data that's based off a Permanent Survey Mark (PSM) (details below). How do I create a CRS so that the data can be brought into GDA94 (4283)?

PSM: 110326
GDA94
Lat: 12D 39M 50.5201"S
Lon: 141D 50M 49.8383"E
E: 591994.116
N: 8599867.976
Zone: 54
Order: 1st Class: A
Adjustment: GSA - WEIPA Control

This seems to be the reverse process of How to define a CRS using a projection rotated about a given point or am I mistaken?

Basically, the question is how to bring a dataset that is fixed to a position that has GPS coordinates (is real world) into real world coordinates as a whole.


texts: Using for-loop to count intersections with bbox

texts:
I'm working on a project that involves counting how many crimes of a certain type happen within 150 meters or 500 meters of a list of addresses. While I've been able to do this successfully without a for-loop, it's a bit cumbersome because of the different crime definitions I need to use, so I'm trying to shorten the process a bit.

The dataframe with address points is about 2000 rows, and the crime points have around 80,000 rows -- definitely a good bit a data, but not so much that it keeps running for 20 mins plus.

Right now, here is the general structure of what I'm trying to do (edited):

function(input dataframe, crimedata1, crimedata2, distance, output){
 for (i in dataframe)
  1. create bounding box row i
  2. calculate number of intersections of boundingbox and file for two different files, based on the year the address is from
  3. return this count
}

countcrimes <- function(data, cdata1, cdata2, d, output){
  names <- c("ID", "year", "geometry")
  output <- data.frame(matrix(nrow = 2071, ncol = 3))
  colnames(output) <- names
  output$ID <- data$ID
  output$year <- data$year
  output$geometry <- data$geometry
  output <- st_as_sf(output, sf_column_name = "geometry", crs = 32616)
  for (i in 1:dim(output)[1]){
      if (output[i,]$year == 1){
      x <- st_intersection(crime1617_v, output)
      output[i,]$n <- st_intersection(cdata1, output[i,]) %>% count()
    } else if (output[i,]$year == 2){
      output[i,]$n <- st_intersection(cdata2, output[i,]) %>% count()
    }
  }
  return(output)
}

vc150m <- countcrimes(bbox150m, crime1516_v, crime1617_v, 150, vc150m)

Right now, I'm dealing with ridiculous runtime and unclear errors -- I'm thinking that the structure of what I'm trying to do isn't working, and I can't quite work out the issue or figure out a better way to structure this process.


texts:
How do I include the joined table when exporting a layer from an MXD using ArcPY.

I am using doing this:

arcpy.FeatureClassToFeatureClass_conversion(layer.dataSource, outLocation,
                                                outFeatureClass)

which does not include the joined information


texts: Combine multiple rows in a CSV using Pandas

texts:
I have over a hundred CSV files that I will join with selections of a shapefile to represent visually in QGIS.
However, in these CSV files, the data is sometimes split into rows of 1A and 1B, or 13A and 13B, for example:

I can't perform the join with rows like this, because in the shapefile the data is never split into 1A and 1B:

For my purposes, I would like to add the content of 1B to 1A (which is different fields of vote numbers) and have a resulting single row (feature) called '1' where the number attributes are now the sums of the respective numbers in 1A and 1B. Another issue with the join of course is that 1A is a 'string', when I want them to be 'int'.
I don't think the 'one-to-many' join works because of the As and Bs.
I'm hoping there is a solution using Pandas that can run through a CSV file and accomplish this merge/summation for every instance of a row split into A and B. Then, I would have this run over a whole folder to perform the same thing on all my CSV files, like this:
for file in os.scandir(src_folder):
    if file.name.endswith('.csv'):
        #perform the merge-row code


texts:
I have the following raster (shp file using readORG):

dput(summary_grid)
structure(list(class = structure("SpatialPolygonsDataFrame", package = "sp"),
    bbox = structure(c(4346000, 3819000, 4445000, 3867000), .Dim = c(2L,
    2L), .Dimnames = list(c("x", "y"), c("min", "max"))), is.projected = TRUE,
    proj4string = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs",
    data = structure(c("Min.   :0  ", "1st Qu.:0  ", "Median :0  ",
    "Mean   :0  ", "3rd Qu.:0  ", "Max.   :0  "), .Dim = c(6L,
    1L), .Dimnames = list(c("", "", "", "", "", ""), "      Id"), class = "table")), .Names = c("class",
"bbox", "is.projected", "proj4string", "data"), class = "summary.Spatial")

as well as the shp file containing lines:

dput(summary_lines)
structure(list(class = structure("SpatialLinesDataFrame", package = "sp"),
    bbox = structure(c(4329488.96922647, 3429159.10800761, 4998503.48859431,
    4055688.10547651), .Dim = c(2L, 2L), .Dimnames = list(c("x",
    "y"), c("min", "max"))), is.projected = TRUE, proj4string = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs",
    data = structure(c("Min.   :0  ", "1st Qu.:0  ", "Median :0  ",
    "Mean   :0  ", "3rd Qu.:0  ", "Max.   :0  ", "Min.   :  448  ",
    "1st Qu.:13229  ", "Median :28235  ", "Mean   :27205  ",
    "3rd Qu.:40724  ", "Max.   :50608  ", "Min.   :  485  ",
    "1st Qu.:13731  ", "Median :29399  ", "Mean   :28635  ",
    "3rd Qu.:43159  ", "Max.   :53607  ", "Min.   :6519754  ",
    "1st Qu.:8741416  ", "Median :9167928  ", "Mean   :8894830  ",
    "3rd Qu.:9414708  ", "Max.   :9762259  ", "Fishing:121  ",
    NA, NA, NA, NA, NA), .Dim = c(6L, 5L), .Dimnames = list(c("",
    "", "", "", "", ""), c("      Id", "   trip_id", "  new_tr_id",
    "     species", "  HELCOM_Gro")), class = "table")), .Names = c("class",
"bbox", "is.projected", "proj4string", "data"), class = "summary.Spatial")

In this shp file with the lines, each lines has a unique identification number called new_tr_id.

My goal is to create a raster file with the number if lines crossing each cell. A line (so the same new_tr_id) can cross the same cell X times, it would be counted X times. I added the figure below as example:

I have some interesting discussions using the package raster but I still can not find a proper answer to this issue.

Some news: I tried the method described by @Spacedman below but it is cutting the lines crossing together:

I don´t have the same result as @Spacedman but I am not that far anymore :)


texts: Arcmap adds vertices to shapefile at dataframe extent?

texts:
I'm running arcmap 10.4.1 trying to make a simple outline from a polygon shapefile. Should be a piece of cake, but a small part of the polygon extends off past the dataframe at the default projection. When I reproject with a 180° central meridian to get the entire thing in one piece, arcmap creates a vertical line where the polygon crossed the dataframe by adding extra vertices.

I can edit these vertices normally such that the line should disappear, but as soon as I deselect the polygon (before I can even save edits) it adds the extra vertices back in!

Has anyone seen anything like this?


texts: Where can I find hydroelectric generation data for the U.S.?

texts:
It is my understanding that the US Army Corps of Engineers has some hydro generation data on their various websites (divisions, districts). For example, there is real-time hydro-electric data generation for the Tulsa District located here: http://www.swt-wc.usace.army.mil/ops_hydropower/hydromainpage.html. However, other districts do not have this same data. Is there any reason why hydroelectric power generation is only reported in some districts but not others? What is a resource I can use to find real-time (and/or even historical) power generation data for other regions (besides the Tulsa area).

I have scoured every single USACE website for all the regions, have examined the various hydro-electric power administrations in the United States (BPA, TVA, SEPA, WEPA) to little avail. Only BPA contained real-time and historical data. I am especially interested in the southeastern part of the United States - states like TN, GA, AL, etc. Proxy generation (like river/lake levels) may also be useful, but I have not been able to find those aside for the Tulsa/BPA areas.

Are there any other resources for this type of data that I am looking for? Even if it is proprietary that would be fine.


texts: QGIS 2.18.10 getting slower after 3 months?

texts:
I'm currently reporting an issue our users are facing in the field.

We are using QGIS / PostgreSQL to collect data on daily basis and since few weeks they are noticing a huge difference in performances with QGIS. For instance: henceforth it takes for them approximatively 10 secondes to open a single form where it took QGIS just half a second to open it at the beginning (3 months ago).
All the QGIS process is slowing down, which is delaying their job and creating issues... they can't afford to wait for the "tablet" in their duty.

They are working on a Microsoft surface pro (i5, 4Go, 128Go SSD) bought in October to replace a previous tablet which was way less performant. At the beginning, the combo surface pro/ QGIS was working extremely well, with high level of performance and reactivity.

The tablet is completely offline, no updates, no connection with internet, no system modifications since the first use. It appears all the others "windows tasks" they need to do (running scripts, opening softwares) are working without any loose of performance so far.

Our QGIS project contains 6 layers, some with geometry, which are only use to open and fill forms in order to fill the database stored in PostgreSQL. The amount of data they have to enter on daily basis didn't change from the time they received the surface pro.

According to web solutions, we :
- unselected all plugins and check : no difference
- selected "show markers only for selected features" and check : no difference

Hence I'm running out of solutions and I'll plan to migrate back to a more stable version like 2.14 and see what it does.

Have you ever experience such a problem and/or would you have some tracks to follow ?

before installing the new version of Qgis (actually I was planning to change it to the 2.16), I tried to create a new project and it worked.

From the feedbacks of the users when we just done it, it seems to be as fast as the beginning for the moment. Nevertheless, we need to wait and see few days before clapping hands..

But why would it have solved the issue ?

any idea/suggestions ?

Simple bug ?


texts: ArcGIS Tools Script Validation...filter list from domain values

texts:
Anyone know if it's possible to build filter list from a set of coded domain values I have in a file geodatabase? I have a script that works taking values from a table and correctly building a filter list but i'm not sure how to access the domain values?

FIELDS = ["Name", "Title", "Organisation", "Address1", "Address2", "Address3", "Town", "PostCode", "Telephone", "EMail"]
TABLE = (r"P:\Photo-Share\New_Projects_Database.gdb\Customer")

rows = arcpy.SearchCursor(TABLE)
self.params[3].filter.list = sorted(list(set(r.getValue('Organisation') for r in rows)))
del rows


texts: How to color a WMS layer using database value

texts:
I have a transparent WMS layer of polygons. It is a map of a city and i want to color the map according to crime rate. Rate values are obtained from database. I dont know how to select a single polygon from that layer and color it according to the rate. Please help. Am trying to implement it using GeoServer and Openlayers and this is a client side operation

This is my map and i want a output like shown below


texts: Automate Parameter Loading (for Labeling)

texts:
I'm doing sidewalk inventory collection for multiple cities throughout the state at times - and the collection process is basically broken down into different criteria when inspecting the walk - so I usually end up with 10 or more different shapefiles (e.g. Poor Condition, Structure, Slopes) per city. With the information gathered, I produce PDFs using Data Driven Pages.

I use Labeling Manager and do Copy Parameters/Paste Parameters from previously loaded shapefiles to populate the expression and symbolize my newly imported shapefiles.

Since the process of moving the files into ArcMap and populating the label expression is always the same, I'm wondering if there is a scripting solution to automate this process?

Initially I thought of ModelBuilder but I can't find any tool that deals with modifying label expression of a feature class - and I'm lacking in the programming department. What I'm envisioning is a menu that open up and allow batch input of feature classes and churns them out with the correct label expressions and symbology loaded.


texts: Better performance for geohash aggregation with PostGIS


texts: Overwrite layer of hosted feature service but not related table in Python

texts: I have a hosted feature service with a related table. The layer in the feature service needs to be updated every month or so. However, the data in the related table needs to remain unchanged. Is there a pythonic way to update the layer without altering the related data? I looked at creating an .sd file with Python but I don't anticipate much success based on what testing I've done. Would deleting the layer's contents and updating through the REST be feasible when I'm dealing with +70 thousand spatial objects?


texts: QGIS not recognising field names in first row

texts:
For some reason, when I add a CSV as a delimited text layer QGIS is not recognizing field names in the first row. This has only been happening the last week or so. I am pulling my hair out.
Below is a screenshot of the csv and a screenshot from QGIS. As you can see QGIS is actually reading row 2 (the first row of data) as the field names.


texts:
I'm currently creating a project for management of highway drains in QGIS which involves the recording the feature attributes and logging of multiple inspections against a feature. I have a layer for the highway_drains and a separate layer for highway_inspections, utilizing the RELATIONS function under project properties works perfectly for the Parent/Child relationship of highway drains to inspections.

I used the drag and drop designer for both layers to reference lookups and style the form. However in order to view/edit the inspections within the form for highway_drains, I have to switch the highway_drains layer from 'drag and drop' designer to auto generate. This results in the lose of the form styling for the highway_drains layer, namely the TABS and order of fields (lookups and field aliases are preserved), styling for highway_inspections is preserved.

My Question: How can I preserve the form style for highway drains?

I'm using QGIS version 2.8.2


texts: Does SAGA in QGIS processing can support no ascii name?

texts:
From Qgis, I am calling through Processing the Saga algs "clipgridwithpolygon"

raster_inside_poly = processing.runalg("saga:clipgridwithpolygon",
                my_raster_temp,
                my_edges,
                my_new_raster)

Since my_new_raster name is all ascii. It's fine.
When no ascii is presented to Processing and SAGA crash with "Encode exception"

My Question :
is this a limitation on the SAGA side or on Processing ?
Someone can propose a solution ?

Problem experienced on Win8 Qgis 2.12.2 Processing 2.12.2 Saga 2.1.2


texts:
I'm trying to select streams that are a certain distance away from cities. I have been able to do it for a few, but I would like to do it to the entire column that contains the cities name without having to do each individual city.

This is what I have so far

CREATE TABLE public.streams_cities AS
SELECT DISTINCT ON (b.gid) b.gid, b.name, b.geom, s.city_name
    FROM public.cities_coutnies s
        LEFT JOIN public.streams b on ST_DWithin(s.geom, b.geom, 10560)
    WHERE s.city_name = 'ARVADA'
    OR s.city_name = 'BRIGHTON'
    OR s.city_name = 'COMMERCE CITY'
    OR s.city_name = 'DENVER'
    OR s.city_name = 'EDGEWATER'
    OR s.city_name = 'FEDERAL HEIGHT';

There must be a less time consuming way.


texts:
I am creating parcel map books for Wisconsin counties. As one might be able to guess, many properties are very large, some spanning multiple sections.  On the other end there are many small parcels in towns and cities.

Using a single sized Index means my parcels are way too small, or I have many pages containing only one property.  Is there a way to create data driven pages that decides which index to use based on parcel count?

My current idea is long and not very efficient both with memory and time, it looks like this:

Create 3 map books, one with Sections, one with Quarters, and one with
Quarter Quarters as the index.
Export to these to PDF's
Back in ArcMap, loop through each of the sections.
Count the number of parcels in the section.
Based on the count, decide whether to grab the Section mapbook page, 4
Quarter map book pages, or 16 Quarter Quarter mapbook pages.
Save these pages in a new folder, when the loop is completed, I'll have my
mapbook.

I am wondering if anyone has had a similar case when working with data driven pages and what their work around was.


texts: What is the gdrive callback url for a self hosted cartodb installation?

texts:
I have a custom install of cartodb running on my own server. I need to integrate Google Drive with it.

I have set up an application on the Google Developers Console, but do not know what to insert for the callback url.


texts: Uploading Layers on ArcGIS Online displays unexpected line?

texts:
I digitized a cyclone track on ArcGIS as a polyline. Then i created a buffer of 100 km using the buffer tool. My parameters when using the buffer are:

Distance- Linear Unit of 100km
Side Type- Full
End Type- Round
Method- Planar
Dissolve Type- None
Coordinate System- Projected (Fiji Map Grid 1976)

Then i wanted to upload it on ArcGIS online. I used the Add layer from file tool. The buffers werent showing so i had to play around with the visibility range of each feature. But when the buffers finally appear, it has a line at the very bottom as shown on the picture below marked blue.

I dont know where that came from.

Any ideas?


texts:
I am joining some tables to a shapefile and exporting with the following code:

import arcpy, pandas as pd, os
import arcpy.mapping
arcpy.env.overwriteOutput = 1
arcpy.env.qualifiedFieldNames = False

arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'F:\Sheyenne\Pixel_Regression\spatial_ts\NDII\monthly_regress'
shapefile=r'F:\Sheyenne\Pixel_Regression\raster_to_point\Sheyenne\NDII\sheyenne_19840517.shp'
out=r'F:\Sheyenne\Pixel_Regression\spatial_ts\NDII\monthly_shapefiles'
if not os.path.isdir(out):
    os.mkdir(out)
#loopo through tables
for table in arcpy.ListTables("*"):
#  #Create a new feature layer, this will remove previous join
  arcpy.MakeFeatureLayer_management(shapefile, "parcelLyr")
#   Join table to feature layer
  arcpy.AddJoin_management("parcelLyr", "POINTID", table, "Month", 'KEEP_COMMON')
  tableName = table.split(".")[0]
  if not os.path.isdir(out):
            os.mkdir(out)
#   #Export joined layer to new feature clas
  arcpy.FeatureClassToFeatureClass_conversion("parcelLyr", out, tableName)

Each table has 290,000 rows (and therefore joins), so maybe there is no faster way to do this.  To join 5 tables with this method takes around 15 hours though which has me wondering if there is an alternative method.

EDIT:

I have explored using cursors to do this.  To start I will show what my data looks like (very simplified).  The shapefile I am joining to (hereafter called targshp) has an attribute table which looks like this:

FID   Shape   POINTID   GRID_CODE
0     Point   1         -999
1     Point   87275     -999
2     Point   87276     -999
3     Point   87278     -999

my csv file (hereafter called joinshp) looks like this:

Month    RGR_Slope   RGR_pvalue   RGR_significant
87275    0.023       0.5979       No
87276   -0.021       0.0139       Yes
87278    0.061       0.2345       No

My join fields are POINTID for targshp and Month for joinshp.  joinshp does not have as many fields as targshp and therefore this would be a type of KEEP_COMMON join.  I am using the following approach using the cursors:

import arcpy
arcpy.env.overwriteOutput = 1
arcpy.env.qualifiedFieldNames = False

targshp=r'F:\Sheyenne\Pixel_Regression\spatial_ts\RGR\cursor_try\sheyenne_19840517.shp'
joinshp=r'F:\Sheyenne\Pixel_Regression\spatial_ts\RGR\monthly_regress\05.csv'
joinfields=['Month', 'RGR_slope', 'RGR_pvalue', 'RGR_significant']
joindict={}
with arcpy.da.SearchCursor(joinshp, joinfields) as rows:
    for row in rows:
        joinval = row[0]
        val1 = row[1]
        joindict[joinval]=[val1]
del row, rows

#
#
arcpy.AddField_management(targshp, "Month", "LONG", "5")
targfields = ['POINTID', 'Month' ]
with arcpy.da.UpdateCursor(targshp, targfields) as recs:
    for rec in recs:
        keyval = rec[0]

        if joindict.has_key(keyval):
            rec[1] = joindict[keyval][0]
        else:
            rec[1] = 0
        recs.updateRow(rec)
del rec, recs

This code runs through without error, and Month is added to targshp but the rest of the data is not actually joined into it, so I am still slightly off here.


texts: How to generate list of polygons that fall within radius circle using Google Maps API?

texts:
I have created a Web app using Google Fusion Tables, the Google Maps API and Derek Eder's Searchable Maps Template that allows someone to enter an address and distance, and then see which US Congressional Districts fall within the radius circle. A blind site visitor asked if it's possible to include a text list of those polygons that intersect with the radius circle. I thought that was a great idea and would like to help visually impaired users make use of the application.

Derek shows how to do a text list of points within the circle, but not intersecting polygons. Would anyone know how to modify this results list code to show intersecting polygons?

I was trying to play around with this code, but unfortunately I'm not very proficient in JavaScript (I use R to prepare all the data).


texts: Naming polylines in sequences from left to right (west to east) using QGIS

texts:
I would like to give each of these polylines an "ID" e.g. 1,2,3,4 etc. from left to right (West to East) or more accurately off a bearing of 125°.
How can I achieve this in QGIS?
Data: https://we.tl/t-SzWkhwRO0D


texts: Applying display properties similar to using ApplySymbologyFromLayer tool in ModelBuilder?

texts:
I have created a layer file using ArcGIS 10.3 ModelBuilder.

I have applied symbology to this layer using the ApplySymbologyFromLayer tool in ModelBuilder.

The tool applied the symbology to my created file, but it did not transfer over the display properties (the options found under the 'Display' tab within Layer properties).

In particular, I wanted the selected resampling method (bilinear) from the source layer (the one whose symbology was copied) to be applied to the ModelBuilder layer file the way the symbology was applied. The model actually creates a different layer for each iterated value, so I'd rather not do it manually for each output.

Is there a way to apply source layer display properties to a newly created layer using ArcGIS 10.3 ModelBuilder?


texts: Unique Text Element per Data Driven Page

texts:
I'm working with data driven pages and need to create a unique text element on every page.

Each text element you add to a data driven page gets added to every subsequent and previous page.

Is there a way around this?

The workaround I've been using is to create the unique note on each page using Bluebeam after the page has been printed which is obviously not ideal.


texts:
I insert the values from lat/long into the points. Those users move around and send me new Lat/long  and of course , I made an update of the position.

I need to find people that are in a distance less than 100Km, which is the best approach to do that?.


texts: How can i create a context menu on Arcgis Desktop.And after creating it,how can i use it while a base command is open and active.I use C# but i can understand visual basic.


texts: European Mountains Shapefile Data Sources

texts: I am trying to locate a GIS freely downloadable shapefile (preferably polygon) of European main mountain chains (not an elevation data), eg Alps, Carpathian Mountains. I tried variuos sources , including Natural Earth, EEA, and many others, but couldn't find the right data. Can anybody help??


texts: Making one field/values equal to another field/values

texts:
I have created a new field called "SFX" in a shapefile called "myshape.shp":

input = r"C:\files\myshape.shp"
arcpy.AddField_management(input, "SFX", "TEXT", field_length=100)

There is now a field called "SFX" but there are no values in the cells. Now, in that same shapefile, I have another field called "SUFFIX".

fields = [field.name for field in arcpy.ListFields(input)]
for field in fields:
    if field == 'SUFFIX':
        print field
>>>'SUFFIX'

Lets say 'SUFFIX' has values 'C', 'D', 'G', etc.... What I want to do, is for each row in 'SFX' I want that cell value to be equal to the corresponding cell value of 'SUFFIX'.

Without making complicated lists using a cursor, how would I do that using UpdateCursor or the Field Calculator?


texts:
I have a database of GPS points. There aren't any tracks, only points. I need to calculate some value for every 100 meters, but sometimes GPS gave a wrong coordinates that lies far from real GPS points, and instead of calculating values for a small square, I have to calculate it for a really big rectangular area.

What is the best algorithm to filter wrong GPS points?

I made a screenshot to help understand:


texts: Is it possible to embed directly a SVG code in an "image element" in the map composer in QGIS?

texts:
I'm trying to implement custom SVG files in the map composer, I was able to add everything I wanted through HTML elements, but I want to add a custom "functional" North Arrow so I need to use the "image element".
PS: I know I can use a URL but I want to avoid using that method.


texts: R-Leaflet: Plotted points on wrong country

texts:
I am learning the leaflet package and am running into a snag. To learn the package I am using the Snow deaths dataset from the package HistData which contains historic Cholera deaths from an epidemic.

When I run the following code the deaths plot correctly put the map layer underlying the dots is Africa rather than London. Any tips on correcting this?

# John Snow Cholera data
library(leaflet)
library(HistData)
library(sp)

data(Snow.deaths)
coordinates(Snow.deaths) <- x+y
proj4string(Snow.deaths) <- CRS("+init=epsg:3109")

leaflet(Snow.deaths) %>% addTiles() %>%
  addCircles(lng = ~y, lat = ~x, weight = 1,
             radius = ~3000, popup = ~case
  )


texts: Pattern fill for raster contours in QGIS

texts: How do I fill the area within raster contours with a pattern fill, instead of colors in QGIS? Is this even possible?


texts: Sharing shapefiles across network for QGIS?

texts:
I am in a very small office of only three computers using QGIS. We wanted to start using one central repository for our datasets as dropbox is used for main network sharing, but would be to slow syncing large datasets.

So I created a simple homegroup which everyone joined at put all the datasets on my free hard drive.

But when I create a project and save it, and then others open it up it is still asking for links to the shapefiles every time.

I realise it is an inelegant solution at this time!

But I just cant work out while the shapefile paths aren't at least sharing across the computers in the homegroup?


texts: Conda Forge or Conda Channel for GDAL on python2.7 and python3.x

texts:
I have two conda environments: One based on python 2.7 and the other based on python 3.5. I want to create a codebase (frequently used functions) that works with both Python versions. There is one module that gives me problems and that is GDAL. GDAL is somewhat more complex because it has a dependency on a C library and is not pure pythonic.

apt-get install python-gdal -y
apt-get install python3-gdal -y

I tried many ways to install GDAL for python and the command line tools and I was wondering which channel is preferred: Conda-forge or the default conda channel? Where to find more information about this and how to make this decision? Mixing channels is not a good idea as far as I understand since GDAL uses a C library.

In addition to the conda install statements I have also seen:

apt-get install python-gdal -y
apt-get install python3-gdal -y
apt-get install gdal-bin

Why are these required?


texts: Zoom (with mouse scroll wheel) into specific area centred by the mouse cursor in MapInfo?

texts:
In MapInfo Professional V11.0, on 64 bit windows 10 OS, when you zoom in using your mouse scroll wheel, it only zooms into the centre of the map shown on the screen. So if you want to zoom into a specific area, you have to zoom a bit, then pan, then zoom some more, then pan a bit more, all to ensure that the area you want to zoom into is in the centre of the screen.

But in other software, like QGIS for example, when you zoom in using the mouse scroll wheel, all you need to do is put the mouse cursor over the area you want to zoom into and then scroll.

How do you achieve this mouse centred scroll capability in MapInfo?

Edit:

As per comments, I tried Map preferences -> display but this is my window and it doesn't include the necessary option.


texts: create grid cells of fixed dimension around a point acting as a centroid

texts:
I've been trying to solve a problem for the last few days but not having much luck.  I'm hoping someone here can point me in the right direction.  I'm new to ArcGIS so I hope my explanation makes sense.

Here's my problem.  I have a layer of points with x,y co-ordinates.  They are laid out in a mecrator projection.  What I need to do is create another layer with a grid.  Each grid cell needs to be 2 degrees high, by 2 degrees wide, with the point acting as a centroid (so the boundaries are at a distance of 1 degree above the point, 1 degree below, 1 degree left and 1 degree right of the point).   The grid cells cannot overlap.  Because of the mecrator projection, I would expect the height of each cell (latitude) to remain the same, but the width (longitude) to expand the further away from the equator they are.

I've tried the thiessen polygon tool, creating buffers, the fishnet, and the minimum bounding geometry tool but its not quite what I'm after.

Any ideas?  Any help is greatly appreciated!

Thanks very much

PS: I've popped a (badly drawn) image up on to G+ of what I want. See here. Figure 1 is what I have, and Figure 2 is what i need to create


texts: Using feature templates for other layers than original to change target layer in ArcMap?

texts:
I was wondering if there is, in ArcMap 10.6, any way to use feature templates in a different layer than the one they where created for.
I have around 10 layers with similiar atribute tables (same columns) and around 30 different possible types of polygons (landuse like forest, meadow, street...). Every feature that I create gets two attributes: landuse name and landuse ID. I created a feature template for each landuse type in order not having to fill in these atributes manually every time. I would just choose the right landuse template and create the polygon.
However i couldnt find a way to use these templates on any other layer (which has the exact same structure). "Copy" only dublicates the template within the layer. When choosing the template from the first layer in the "create features" dialog, the new feature is also written to that layer...

In other words: is there a way to change the target layer of a feature template?


texts: Using local background map in ArcGIS Mobile?

texts:
My customer requirement is to support a MrSID image directly in ArcGIS Mobile.

The details are sketchy but I think the intent is that they want a local background raster Map while the ESRI.ArcGIS.Mobile.Map gets used to handle other tasks (?vector display, updating GDB?).

The Map control seems designed to consume Map Services via the MapCache object.  Is there any way to get the MrSID image to render, pan and zoom in ArcGIS Mobile?

A "How to" reference is perhaps too much to hope for.  But I'll appreciate any related suggestions. For example, how might we use a local TIFF or JPEG image as the background map?

The control that renders for us in ArcGIS Mobile is ESRI.ArcGIS.Mobile.Map.  This control exposes a “Layers” property which is a collection of MapLayers. If we make a custom MapLayer that knows how to decode a MrSID file, it can read a local file, decode it and blit the pixels to the screen via it's overriden Draw() method.

The difficulty, of course, is that there is currently no Native ARM or pure Managed SDK that will do this.  But if we work through that, then we could solve the integration problem.  The 9.3.1 Mobile SDK includes a sample called LiveTraffic that illustrates the technique:  derive a custom MapLayer and override the Draw() method().


texts: Selecting features by location based on selected feature using ArcObjects?

texts:
I am very new to ArcGIS and ArcObjects.

I am writing a code for a button that will detect a selected feature and then select other features based on the feature that I had selected please take a look at my code

public class LandmarkBtn : ESRI.ArcGIS.Desktop.AddIns.Button {

    public LandmarkBtn() {
    }

    protected override void OnClick() {
         var doc = ArcMap.Document;
        var focusMap = doc.FocusMap;
        IEnumLayer layer = focusMap.get_Layers(null,true);
        ILayer lyr = null;
        IFeatureClass Lmarks = null;

        while ((lyr = layer.Next()) != null) {

            var featureLayer = lyr as IFeatureLayer;
            if (featureLayer == null)
                continue;

        var name = featureLayer.FeatureClass.AliasName;

            if ("landmarks".Equals(name.ToLowerInvariant())) {

                Lmarks = featureLayer.FeatureClass;
                break;
            }
        }

        if (Lmarks == null) {
            MessageBox.Show("Sorry, No LM Layer Found");
            return;
        }

        var selection = (IEnumFeature)focusMap.FeatureSelection;
        IFeature ftr = null;
        ftr = selection.Next();

ftr now holds the selected landmark(feature), I want to perform a query that will select features from another layer that are within 100 meters of the landmark feature, can you please direct me onto where I should look or provide me with a sample code that might help?


texts:
I'm trying to create new vector layer based on existing. Suddenly long  attribute names are truncated. Here is my code:

@staticmethod
def CreateNewEmptyLayer(sourceLayer, newLayerFileName):
    provider = sourceLayer.dataProvider()
    writer = QgsVectorFileWriter(newLayerFileName,provider.encoding(), provider.fields(), sourceLayer.wkbType(), provider.crs(),"ESRI Shapefile")
    if writer.hasError() != QgsVectorFileWriter.NoError:
        print("Error when creating shapefile: ", writer.errorMessage())
        return False;

    del writer;
    return True;

Source layer is MapInfo TAB file.

Added code to print info for source layers:

for field in provider.fields():
        print field.name(), field.typeName()

and result is ->source layer atrtribute names are NOT truncated to 11 symbols. So interesting is why QgsVectorFileWriter is truncating them?


texts: Creating DEM from topography map in QGIS

texts: For my geology diploma project I need to create  DEM from topography map my professor gave me to use as a base map. The map is georeferenced. My idea was to draw lines over the izolines of the map, then specify thair elevation and then create DEM. Equidistance of the izolines is 10m. I don't know how to do it or is it even possible in QGIS.


texts: How to print ArcGIS tile map using MapFish print


texts: Running gdal_contrast_stretch command in Windows

texts:
I am trying to run a GDAL command that converts the bit depth while applying a stretch function (gdal_contrast_stretch -ndv 0 -linear-stretch 70 30 ./tmp/b2-projected.tif ./tmp/b2-8bit.tif). However, I can't find a way to do it in MS Windows. I have checked the precompiled GDAL binaries in GISinternals and they don't carry the above mentioned command.
If not available, what do you guys suggest to use to get the same output using CLI tools in Windows?


texts: QField unable to load compressed TIFF files

texts:
I have come across the following error in QField when trying to load transformed and compressed GeoTIFFs. All uncompressed TIFFs seem to load fine.

  The following layers could not be loaded, please review and reconfigure the QGIS project

All the paths in my project are relative (as recommended) so I am not sure why I am receiving this error. The raster data is compressed via the GUI (tested with asc, xyz and TIFF to TIFF):

Right Click -> Export -> Save As -> GeoTiff with a high compression

I understand the official recommendation (https://qfield.org/docs/project-management/dataformat.html#raster-data) is to use a geopackage raster with pyramids, however my goal is to create a workflow for people with limited GIS understanding/patience (forestry). So far compressed TIFFs are easy to create and function well in QGIS.


texts: How do you hide a feature within a point layer?

texts:
I know this is basic, but for the life of me I can not remember how to hide/disable a specific feature point within a point layer. This seems like it should be a 2-3 click solution, but working with Arcmap I should always know better. One would think ESRI could simply make a right-click/hide point (polygon, or line), then have a "show all points in layer option" to bring it back.
Thanks for the help!


texts: Table for feature layer in PostGIS database cannot be published and enabled on GeoServer

texts:
I have tried to publish a table for feature layer in a PostGIS database to GeoServer 2.2.5. The layer was added and published. However, the layer has got a yellow exclamation mark in the "Enabled?" field in Layers view. It is also missing from the Layer Preview page.

Everything seemed OK during the publishing process. I could also add the table to QGIS and it's shown well. Anything I was doing wrong?

Can anyone give me some hints about it? Thanks!

Alex


texts: Trouble "extracting by mask" with a DEM layer and a shapefile?

texts:
I'm using ArcMap 10.2 and I'm currently trying to clip a county-wide DEM layer (.lyr file) down to a smaller area using a shapefile. The data frame coordinates system is set to US Feet which is what I need. The DEM's coordinate system is in meters and I cannot change that. The shapefile I'm using to clip the DEM I have tried both setting to Meter and to US Feet. I'm using the Extract by Mask tool in the Spatial Analyst > Extraction toolbox. Each time I receive the same error message: Output Raster is listed as "empty" (I do have a workspace set up locally) and the error codes I receive are: 010429, 010213, and 010067. I've looked these up but I know I've gotten this process to work before, I think it's more a coordinate system problem than anything.

Here is the error message I've received each time after trying to extract:

I've found and downloaded a raster/DEM that is in US Feet. I've set the data frame to US Feet, and I've made sure the shapefile to mask is set in US Feet. In the "Extract by Mask" tool window, I went into Environemnts and set the Output Coordinates to match the DEM and the Raster Analysis Cell size to match the DEM. The workspace is local and I've tried decreasing the size of the shapefile I'm using to extract (incase it was just too big). I'm still getting the same three errors. From what I can understand, the first error (010429) seems to indicate that arcmap can't read the DEM layer correctly? It doesn't recognize it somehow? I've downloaded and used DEM's from three different providers so I don't think it's a sourcing issue.


texts: Total count of values at each pixel in a raster stack/ 3d array

texts:
If you have an array with x, y and z dimensions, how can I create an array of x, y dimensions that holds a count of how many z values at each location (x,y) are noData?

I would prefer to use numpy and python for this task, but could also use arcpy (I am working with raster data) or R.


texts: How can I prevent markers from being clipped?

texts:
I have a series of point that I want to draw circles with varying radii around. The radii depend on the zoom level and represent a 1.1 mile circle from the point. I'm running into an issue where the circles will draw but are clipped.

All I've been able to find so far is to increase the buffer size on the source. I've maxed this out at 64 and it doesn't make a difference.

Is there any way to be able to draw markers around points with this kind of large size? Some other approach that I'm missing or another buffer I could adjust?


texts: I have two layers in a LayerGroup that I can toggle between and would like to update a control(bottom left in example) with the grades of the currently selected layer. How can I go about doing this in leaflet? This is what I currently have here


texts: Statistics of NDVI mean by shapefile using GEE

texts:
I am new researcher using GEE and want to use a shapefile (e.g.,Provincial boundaries) to achieve NDVI mean of each province using GEE, how should we do it?

code:

// Create a mosaic of Landsat 5 images from June through September, 2007.
var allBandMosaic = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .filterDate('2016-03-01', '2016-11-30')
  .select('B[0-7]')
  .median();

// Var geometry = ee.Geometry.Rectangle([77, 17, 135, 56]);

// Use the normalizedDifference(A, B) to compute (A - B) / (A + B)
var ndvi = allBandMosaic.normalizedDifference(['B5', 'B4']);

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

// Center the map
Map.setCenter(103.16, 34.82, 8);

// Display the input image and the NDVI derived from it.
Map.addLayer(ndvi, {min: 0, max: 1, palette: palette}, 'NDVI');

// Export the image, specifying scale and region.
//Export.image.toDrive({
//  image: ndvi,
//  description: 'ndviToDriveExample',
//  scale: 30,
//  region: geometry
//});
var fc = ee.FeatureCollection('users/qiangsun/shengjie')
  .filter(ee.Filter.stringContains('ENAME', 'Gansu'));
Map.addLayer(fc,{color:'FF0000'},'shenjie');

var means = ndvi.reduceRegions({
  collection: fc,
  reducer: ee.Reducer.mean(),
  scale: 30,
});
means = means.get('ndvi')
print('gansu',means);

However, I want to batch process for all provincial mean NDVI in China, and why the result for Gansu is null？.


texts: PostGIS ST_Intersection of polygons can return lines

texts:
When clipping the polygons of one table with polygons in another, ST_Intersection can return a set of results that can be handled with ST_Dump. The returned multiple geometries are not necessarily ST_Polygon but also ST_LineString (probably also a point). So when running a query

INSERT INTO c (geom)
  (SELECT (ST_Dump(ST_Intersection(a.geom,b.geom))).geom
  FROM a INNER JOIN b ON ST_Intersects(a.geom, b.geom));

when trying to populate the table "c" with clipped polygons, it fails with
ERROR:  Geometry type (LineString) does not match column type (Polygon)

I did another nested SELECT statement so only the polygon geometries came through, like:

INSERT INTO c (geom)
  (SELECT geom FROM
    (SELECT (ST_Dump(ST_Intersection(a.geom,b.geom))).geom
    FROM a INNER JOIN b ON ST_Intersects(a.geom, b.geom))) AS cl
    WHERE ST_GeometryType(cl.geom)='ST_Polygon');

As this is a bit cumbersome I am wondering if there is a more elegant solution for dropping invalid geometries?


texts: Donut Polygon remove in Building SHP exported from Autocad

texts:
I have created Buildings footprint in Autocad and exported to SHP which is having so many Donut polygons. Is there any way to clip polygons in SHP files (Batch)

I have found one script on net but it is creating issue as it is clipping overlapping errors too.

Modified Script:
import arcpy
import pythonaddins
import os ,sys

System Variables:

mxd = arcpy.mapping.MapDocument(r"Current")
df = arcpy.mapping.ListDataFrames(mxd)[0]
arcpy.env.overwriteOutput = True

Local Variables:

Vendor = arcpy.GetParameterAsText(0)
donutLayerUnion = arcpy.GetParameterAsText(1)
idField = arcpy.GetParameterAsText(2)
Donut_Processed = arcpy.GetParameterAsText(3)

Union analysis

arcpy.Union_analysis ("Vendor","donutLayerUnion","ALL","","GAPS")

Delete Identical analysis

arcpy.DeleteIdentical_management ("donutLayerUnion","Shape;Shape_Area","","")

Copy Feature Class

arcpy.FeatureClassToFeatureClass_conversion ("donutLayerUnion","C:/Temp","Donut_Processed.shp",where_clause="",field_mapping="",config_keyword="")

Delete Layers

arcpy.Delete_management("c:\temp\Donut_Processed.shp", data_type="ShapeFile")

Delete layer from active frame

for df in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "donutLayerUnion", df):
        arcpy.mapping.RemoveLayer(df, lyr)
pythonaddins.MessageBox ('Process Done..','Kailash',0)


texts: Selecting parts of a polyline on a map leaflet

texts:
Is it possible somehow to select parts of a polyline on a leaflet map?
Say, I have a star-formed polyline, and I'd like to select one peak of it (and then remove it from the polyline).


texts:
The following script is supposed to find and replace all shapefiles in a mxd (that exist in the folder) to the new SDE data source. No errors when running the script however, when I open the MXD in ArcGIS 10.1SP1 nothing changed.  Any suggestions what is wrong?

import arcpy, os, glob

#find all the MXD's in the directory tree
for root, subFolders, files in os.walk(r"C:\Users\lmuhammad\Documents\ArcGIS\SearchandReplace\Testing"):
for filename in files:
    fullpath = os.path.join(root, filename)
    basename, extension = os.path.splitext(fullpath)
    if extension.lower() == ".mxd":
        print "------------------------------"
        print filename
        #open the map document
        MXD = arcpy.mapping.MapDocument(fullpath)
        #get all the layers
        for lyr in arcpy.mapping.ListLayers(MXD):
            #get the source from the layer
            if lyr.supports("datasource"):
                source = r"R:\Shapefiles\Road_Casings\RoadCasing.shp"
                print "%s -> %s" % (lyr, source)
                basename, extension = os.path.splitext(source)
                if extension.lower() == ".sde":
                    #This is the NEW SOURCE that you want to point to
                    datapath = r"Database Connections\Connection to GISD.sde\GISDEV.TRANSPORTATION\GISDEV.PAVEMENTEDGE"
                    #replace the old path wih the new
                    lyr.findAndReplaceDatasource(source, datapath, "NEWROAD", False)
        #save your changes
        MXD.saveACopy(r"C:\Users\lmuhammad\Documents\ArcGIS\SearchandReplace\Testing\test2.mxd")
        del MXD


texts:
I am trying to use the WMS and WMTS services of the Civil Protection Department weather radar:
https://dpc-radar.readthedocs.io/it/latest/services.html
However, in QGIS 3.18 it is not known that the images are not updated when they should do it every 5 minutes as per web viewer:
http: //www.protezionecivile.gov.it/radar-dpc /
With F5 nothing happens and I also tried not to set the use of caching in display.
I tried with this code without success
from qgis.core import QgsProject
layer = QgsProject.instance().mapLayersByName('vmi')[0]
layer.dataProvider().reloadData()
layer.triggerRepaint()


texts:
I've got a working example of a map that adds a marker layer with the creatLayer method. I am trying to add a click listener to the carto markers. I have tried the examples here and here to no avail.

Here is where I am at right now. You'll see the markers are loading, but no interaction on click (it should call an alert). I have tried a lot of variations on the click listener code at this point.

Code pasted below as well.

<!DOCTYPE html>

<title>How to add markers to a CARTO map</title>
<meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
<meta http-equiv="content-type" content="text/html; charset=UTF-8"/>
<link rel="shortcut icon" href="http://cartodb.com/assets/favicon.ico" />

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.15/themes/css/cartodb.css" />


texts: How can I create map tiles for an irregular shaped area, like a state?

texts:
I'm using Ubuntu 14.04, OSM in a Postgres database, and Tilemill to style and export to a mbtiles files. How can I generate tiles for an irregular-shaped (not rectangular) area. For example, if I want to generate tiles for the US states of Michigan or Illinois, I currently specify a rectangular bounding box based on the min and max latitude and longitude. This results in an enormous amount of unnecessary tiles being generated for, in the case of Michigan, Wisconsin and Canada.

If it can't be accomplished with my current tools, can I alter my process to incorporate a tool that allows me to use an irregular area?

Is there a tool I can use to trim the source OSM data to an irregular boundary? Even if I'm limited to a rectangular box when generating tiles, at least there will be no data and hopefully I get blank tiles or a single tile repeated over the areas with no data.


texts:
I have been given an .slpk file generated by Drone2Map. I can import the .slpk file into ArcGIS Pro but I'm unable to edit the 3D layers in the package.

I have not been able to find a way to extract the layer data from the .slpk file so that I can make it editable. Is there a way to do this?


texts: Transforming map with "custom" CRS to WGS84?

texts:
My map opens with QGIS saying "This layer seems to have no projection specification" and generates the following CRS:
+proj=aea +lat_1=31.5 +lat_2=24 +lat_0=24 +lon_0=-84 +x_0=400000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Looks like some form of Lat/Lon coordinates.  How can I transform to WGS?


texts: Raster image does not show up in Data View

texts:
I'm a new GIS Student. I created and imported 3 raster files for last week's lab. In our lab this week, we are creating a vector image using the raster data. When I import the 3 raster files, it appears to work with no import error. I cannot see any of the 3 but they are in a layer and in the catalog.

ArcMap 10.2

ncols 15
nrows 15
xllcorner 2105000
yllcorner 660000
cellsize 5000
nodata_value -9999
0 0 0 1 1 1 2 0 1 0 0 1 1 1 2
2 0 0 0 0 1 2 0 1 0 0 1 1 2 2
2 0 0 0 0 0 1 1 1 0 0 1 2 2 0
2 0 0 0 0 0 2 1 1 0 1 1 2 0 0
2 2 0 0 0 0 2 0 1 1 1 1 1 0 0
0 2 1 1 1 1 1 1 1 1 1 1 1 0 0
1 1 0 0 0 0 0 0 1 0 2 0 1 1 1
0 1 0 0 0 0 0 1 0 0 2 0 0 0 0
1 2 2 2 2 2 2 1 2 2 2 2 2 2 2
1 0 0 1 1 1 1 1 0 0 0 0 0 0 2
1 2 2 2 2 2 2 1 1 0 0 0 0 0 2
2 2 0 0 0 0 0 0 1 1 0 0 0 2 2
2 0 0 2 0 0 0 0 2 0 0 0 0 2 2
2 2 2 2 0 0 0 0 2 0 0 1 0 2 1
0 1 1 1 1 1 1 1 2 2 2 1 1 1 1


texts: Error mapping data frame in R: is.data.frame(map) is not TRUE

texts:
I am trying to generate a map using the following data and code:

fileURL <- c("https://drive.google.com/open?id=0B-jAX5hT2D3hNnVtLVhROENKRGs")
require(data.table)
require(tidyr)
require(dplyr)
ge.planted <- fread(fileURL, na.strings = "NA")
ge.planted[,-1]
colnames(ge.planted) <- c("region", "type", "crop", "2000", "2001", "2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012", "2013", "2014", "2015")
ge.planted.tidy <- ge.planted%>%
        gather(year, percent.planted, -c(region, type, crop))
ge.planted.tidy$region <- tolower(ge.planted.tidy$region)

require(ggplot2)
require(maps)
states <- map_data("state")
snames <- data.frame(region=tolower(state.name), long=state.center$x, lat=state.center$y)
ge.planted.m <- merge(snames, ge.planted.tidy, by="region", all=TRUE)
ge.planted.m$region <- as.character(ge.planted.m$region)
qplot(long, lat, data=ge.planted.m, geom="map", map = region, fill=percent.planted)

When I run this the following error is returned:

Error in eval(expr, envir, enclos) : object 'region' not found

I initially thought that I was getting the error because region was a factor and not a character variable, but I have converted it to a character variable (second to last line of code) but the error remains the same with or without it.


texts: unable to save qgis as default project or using save as option

texts:
I trying to save qgis default template but i got an error

  Error save current project as default

I am usuing the option of qgis

Setting -> Option -> General-> Set current project as default

please help me why I am unable to set project as default

even When I tried to save as project that also I cant save


texts: Extracting multiple values out of table in ArcGIS ModelBuilder?

texts:
I run the OLS (Ordinary Least Squares) tool in my model and output a coefficient table. Now I need to use the slope and intercept values in the coefficient table in a Raster Calculator. I've tried using the different iterators (Row Selection, Field Values) but they only output one value at a time.

Does anyone know how to get "SearchCursor"-like functionality in ModelBuilder?

Here's a picture of the model:

and a picture of "coef_table":

How do I get the two values in the Coef column as separate values that I can use in the Raster Calculator?


texts: Creating color tables for floating point raster maps in GDAL?

texts:
I've got a script that dumps a floating point raster that I'd like to automatically color.  Currently I load this raster in ArcGIS Desktop or QGIS and manually create a color table every time I re-run my model.  I'd like my model to automatically create a nice color table, but the only functionality I can find is to create a color table for integer only rasters.

My approach so far is to output 2 rasters, one with floating point values rendered in grayscale and a second integer raster that has the colors on it.  I can flip between the two to see values, but it's far from optimal.

Is there a technique to procedurally create floating point color tables in GDAL that automatically load with QGIS or ArcGIS Desktop?


texts: QGIS – Edits to shapefile attribute table won't save special characters

texts:
I'm working in QGIS with a Government of Canada shapefile of all 338 federal electoral districts.

The encoding is UTF-8.

In the attribute table, all of the em dashes ("—") that should appear in the riding names have been replaced by question marks.

I've used the field calculator to replace the question marks with the proper em dashes, but those changes won't save properly. If I close and re-open the attribute table, the question marks re-appear.

But if I replace the question marks with non-special characters (e.g. "AAA"), there doesn't seem to be an issue. So the problem seems to be tied to the em dash.

Any ideas on what might be going on?


texts: How to recode and recalculate Field Values using ArcPy?

texts:
I am trying to figure out how to do several calculations through python and the arcpy set.  I have an large dataset (7+ million records of individuals), and I am trying to code a value of 1 that denotes if the individuals are married.

Now I do not have marriage status, but I might be able to get it through following through this logic...

I have last name and address, I can create a new field that has this new identifier and determine if more than one person with the same last name lives at that address.
I have Age, so I could subtract the highest value from the second highest value to get the difference, and if that difference is say less than 15 years, assume they are married.
I do not have gender, but I do have first name.  I have a large dataset of first names to determine gender, but there are too many names that are for both men and women, and I want to stay away from this as this process will be too involved.

So I can create the new field easy through update cursor.  I dissolve it and join it back to the data set to get the number of people with eh same last name at an address.  But after this, I am running into the my issues.

I would like to subtract the largest number that has the same id field, from the second highest number.  But I am not sure how to write/code that in python to run it in ArcGIS.  Is this even the best solution to solving this problem?  I am welcome to all suggestions to help me get through this problem!!!!


texts:
Following the Arcmap 10.6.1 help content instructions, I'm attempting to create a tile package on my local drive, using File > Share As > Tile Package...

However, I get this error:

I find “invalid geometry” odd, because all of the data in the TOC (seven vector layers) is SDE – and all SDE vector data supposedly has its geometry checked upon ingestion by default.  But just in case, I ran the Check Geometry tool on all my layers, and as expected, there were no errors returned.

All layers are in the same CRS.

Customize > ArcMap Options > Sharing > Enable ArcGIS Runtime Tools is checked on.

The Geoprocessing Results window is empty, so I can’t examine the error further.

When I do an internet search on the phrase The junction has invalid geometry, nothing matching is returned.

Any thoughts?


texts: Calculating forest connectivity for raster

texts:
I'm trying to calculate some sort of forest connectivity value for each pixel in my land cover classification. I would like to do this by calculating how many other forest pixels are directly connected (no gaps) to each forest pixel in the data. Is there a tool or approach for this?

What other means of calculating forest connectivity you would suggest?


texts: Building a rectangle around a font marker using the geometry generator

texts:
In able to show labels in a legend
(see my other question regarding the legend here: QGIS styling a legend)
I need to show points as rectangles with the same height but a length based on the length of text in a font marker.

This is a screen view of the points with the font markers with text of differing length.

I tried the buffer function but have so far only busted my brain because I could not work out how to do it.


texts: Reclass all values < 0 to -9999 in PostGIS

texts:
I have a tiled raster in a PostGIS table and I would like to change the value of all pixels which have values less than 0 to -9999. How should I do this?

I suspect it involves ST_SetValue and possibly a SQL sub-statement, but I'm not sure how to put these together.


texts: Create transparent background for svg marker in QGIS?

texts:
I have an image that I want to use as a SVG marker in Print Composer.

I created it as a .png file in Inkscape and then saved it as .svg. But the background transparency does not come across from the png to the svg.

So I guess I have two questions. Can I adjust the transparency of the svg background in QGIS when I set the icon style? (and by default, how?)

And, if not, does anyone know how I can preserve the transparency when I create the svg in Inkscape?


texts: QGIS Plugin: display dialog for entering username and password to connect to the database

texts:
I do not want the password and username to be saved in the QGIS project.

Instead of standard QGis dialog(QgsCredentialDialog) I want to display "my dialog" for entering username and password.
On the "button connect" :connect to DB(Postgresql) (just like standard dialog of qgis is doing but I want to do it with my dialog)

The idea is to connect the signal that is emitted when the project starts  to load into QGIS...

Anyone have solution for this?


texts: Displaying several data rows in popup using Leaflet for R?

texts:
While trying to map some quantitative data with leafletin R, I manage to place pie-charts with leaflet.minichartson the map, like in the following example:

### data ###
d <- structure(list(
   area_name = c("Alscheid", "Alzingen", "Angelsberg", "Angelsberg", "Angelsberg", "Arsdorf","Asselborn", "Asselborn", "Baastenduerf", "Bartringen"),
   surveyID1 = c(1510L, 24L, 382L, 1429L, 1061L, 526L, 1524L, 1281L, 2169L, 1292L),
   Alter = c("25 bis 34", "25 bis 34", "15 bis 24", "15 bis 24", "25 bis 34", "15 bis 24", "35 bis 44", "35 bis 44", "25 bis 34", "25 bis 34"),
   latitude = c(49.970094, 49.56885, 49.76374, 49.76374, 49.76374, 49.860547, 50.09624, 50.09624, 49.891134, 49.6038), longitude = c(6.007703, 6.16394, 6.1592, 6.1592, 6.1592, 5.842195, 5.97425, 5.97425, 6.164145, 6.0782),
   count_all_variants = c(1, 1, 3, 3, 3, 1, 2, 2, 1, 1),
   var1 = c(1L, 1L, 2L, 3L, 3L, 1L, NA, NA, NA, 1L),
   var2 = c(NA, NA, 1L, NA, NA, NA, NA, NA, 1L, NA),
   var3 = c(NA, NA, NA, NA, NA, NA, 2L, 2L, NA, NA)),
   row.names = c(NA, 10L), class = "data.frame")

### mapping ###
library(leaflet)
library(leaflet.minicharts)

tilesURL <- "http://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}"

title <- tags$div(HTML('
Fussball
'))
basemap <- leaflet(options = leafletOptions(zoomControl = FALSE, minZoom = 9, maxZoom = 10, dragging = T)) %>%
   addTiles(tilesURL) %>%
   fitBounds(6.1, 49.4426671413, 6.1, 50.1280516628) %>%
   addControl(title, position = "topleft")

colors <- c('#7fc97f','#beaed4','#fdc086')

basemap %>%
   addMinicharts(
   d$longitude, d$latitude,
   type = "pie",
   chartdata = d[, c("var1", "var2", "var3")],
   colorPalette = colors,
   popup=popupArgs(
   labels=c("Fussball", "Futtball", "Foussball"),
   html=paste0("
", d$area_name, "
",
               "ID: ", d$surveyID1, "
",
               "Alter: ", d$Alter
               )
   ),
   width = 60 * sqrt(d$count_all_variants) / sqrt(max(d$count_all_variants)), transitionTime = 0
)

Working example

In addition to that, I also would like to populate the popup window for each piechart with all the data for the respective location. For the moment, only one data row is displayed, showing 'area_name', 'surveyID1' and 'Alter'. Taking the location 'Angelsberg' as an example (in the middle of the map), I want the popup showing the data for all (3) data rows, which make up the piechart, i.e.:

Angelsberg

ID: 382 Alter: 15 bis 24

ID: 1061 Alter: 15 bis 24

ID: 526 Alter: 25 bis 34

I assume that I have to pass some kind of list/array to html but I have no clue how to achieve this here.


texts: GIS data - Geoprocessing Service to receive data

texts:
I would like to move 8 shapefiles from one intranet to another intranet via arcgisserver (or an other method).  Creating a FTP site is not an option.  I think I saw somewhere on how to send data to a geoprocessing service.

Is this feasible?  Sending a shapefile to a geoprocessing service outside of the api with say python?  Thanks.

My basic need is to move 8 shapefiles from one intranet to another intranet.


texts: How to get WYSIWYG of raster layer in print composer in QGIS?

texts:
I currently have a base raster layer from OSM Standard used in my project. In the regular editing platform for QGIS, the place identifiers for cities are easily readable.

However, once I enter layout mode to prepare the map, the names become small and unreadable. I think this has something to do with the resampling but I'm not certain how to adjust it.

How can I modify the raster layer in layout mode so the place identifiers can be read, like in the editing mode?

I'm reposting this because the previous question I was directed to had no definitive answer. The recommendation to change the dpi to 96 severely limits the height/width attributes of the exported map as an image. Also, changing part of the plug-in code for OSM did not work. In addition the recommendation to download another program and to modify an XML file in the new program's repository is overly convoluted.


texts: geoserver problem to join sql-server table to a layer

texts:
i have a layer published on geoserver based on a connection to a postgres/postgis Table and i need to add another table existing on sql server to do a data jointure table to table in order to add external data to my layer..

is there any way to do it?

tnks in advance to your Help.


texts: I am fairly new to GIS and QGIS, so bear with me. I have a shapefile with various point features and I need to make a heatmap based on an attribute (i.e. # of individuals per sighting; point feature = 1 sighting). Is there a way to base a heat map on an attribute other than by styling as heatmap?


texts: How to perform thinning/skeletonization morphology operation on GEE?

texts:
In google earth engine (GEE) they provide a morphological operator for erode and dilate as follows (Ref: here)

image
    // Perform erosion
    .focal_min({kernel: kernel, iterations: 2})
    // Perform dilation
    .focal_max({kernel: kernel, iterations: 2})

As of right now I am using ee.Algorithms.CannyEdgeDetector to perform edge detection and aggregate over several channels using mean operation. Then I use OTSU's method for binarization (Ref here) to get final binary edge

However, binary edge is too thick, and I want to perform thinning/skeletonization operation over final edge operation. I did some googling and I am not sure of its implementation on GEE. I found that there exists several algorithms (Ref here, here, and here)

This lead me to my original question on how to perform such operation?

For current example implementation can be found here


texts: raster2pgsql: Strange error when creating a table

texts:
When I try to load an ASC file into my psql setup, I get a strange error message. I have created the SQL input file via

raster2pgsql -s 4236 -I -C -M gis_data.asc -F public.lightning > lightning.sql

The parameters follow exactly this link here. When I try to load the data into pgsql with

psql -d database_name -f lightning.sql

I get the following error message:

BEGIN
psql:lightning.sql:2: NOTICE:  CREATE TABLE will create implicit sequence "lightning_rid_seq" for serial column "lightning.rid"
psql:lightning.sql:2: ERROR:  type "raster" does not exist
LINE 1: ...lic"."lightning" ("rid" serial PRIMARY KEY,"rast" raster,"fi...
                                                             ^
psql:lightning.sql:3: ERROR:  current transaction is aborted, commands ignored until end of transaction block
psql:lightning.sql:4: ERROR:  current transaction is aborted, commands ignored until end of transaction block
psql:lightning.sql:5: ERROR:  current transaction is aborted, commands ignored until end of transaction block
psql:lightning.sql:6: ERROR:  current transaction is aborted, commands ignored until end of transaction block
ROLLBACK
psql:lightning.sql:8: ERROR:  relation "public.lightning" does not exist

I thought that the SQL input file would instruct pgsql to create the required tables, especially since the lightning.sql file starts with

BEGIN;
CREATE TABLE "public"."lightning" ("rid" serial PRIMARY KEY,"rast" raster,"filename"
text);

Do I need to create the tables manually?
What is the required type raster?
Where is the mistake I do to get the GIS data into the pgsql database?

PS: Interestingly, my lightning.sql file contains a strange bit sequence. Is that normal? I thought it would be a full ASCII file with SQL instructions?

BEGIN;
CREATE TABLE "public"."lightning" ("rid" serial PRIMARY KEY,"rast" raster,"filename"
text);

INSERT INTO "public"."lightning" ("rast","filename") VALUES
('01000001001975ADBD4F55B53F1975ADBD4F55B5BF0000000000003BC0B4CA4C69FDFF514000000000000 ... ')
CREATE INDEX "lightning_rast_gist" ON "public"."lightning" USING gist
(st_convexhull("rast"));
ANALYZE "public"."lightning";
SELECT AddRasterConstraints('public','lightning','rast',TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,,
FALSE,TRUE,TRUE,TRUE,TRUE,TRUE);
END;
VACUUM ANALYZE "public"."lightning";

Was maybe the conversion not OK?


texts: Using accumulate function for temperature units?

texts:
I am working on a piece of code that is designed to aggregate growing degree days for a particular region.

To start, the equation for growing degree days is: Daily GDD (°F) = Daily Average Temperature °F - 50 °F

So I am trying to compile a function for MODIS data where:

The code calculates the average temp for each day
Subtracts 50 °F from the average temp for each day
Sums up the total values of the output for a given time range
Places that value on the map

The goal being that when looking at different areas throughout a given time range on a map, you will see variations in the total heat units accumulated.

Here is my code

var dataset = ee.ImageCollection('MODIS/006/MYD11A2')
                  .filter(ee.Filter.date('2015-04-01', '2015-05-20'));
var landSurfaceTemperature = dataset.select('LST_Day_1km');
var landSurfaceTemperatureVis = {
  min: 13250.0,
  max: 15690.0,
  opacity: .46,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};
Map.setCenter(-85.60371794450282,44.73590436363271, 10);
Map.addLayer(
    landSurfaceTemperature, landSurfaceTemperatureVis,
    'Land Surface Temperature');
// Create a geometry representing an export region.
var features = ee.Geometry.Rectangle([-85.1417893413635, 45.31413490213395, -86.125065708551, 44.65070625463291]);

// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: landSurfaceTemperatureVis,
  description: 'imageToCOGeoTiffExample',
  scale: 30,
  region: features,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});

I am completely lost as to A. how to average & subtract x from each day within given time range, and B. how to accumulate those values onto a map.


texts: Writing Non Spatial Table into Geopackage - Python

texts:
Just wondering how to write a non spatial (aspatial) table to geopackage using python. I have attempted using Fiona,

geopackage = 'some_geo.gpkg'
layer = 'test_layer'
schema = {'geometry': {'type': 'aspatial'}, 'properties': [('filename', 'str')]}

with fiona.open(geopackage, 'w', layer=layer, driver='GPKG', schema=schema) as db:
    record = {'geometry': {'type': 'aspatial'},
              'properties': {'filename': filename[:-4]}}
    db.write(record)

This comes up with a Attribute Error because its trying to strip a dict, but I dont think thats the real issue here.


texts: Hardware requirements for one instance Geoserver/PostGIS

texts:
What hardware minimun requirements do you recomend for one instance of geoserver and postgis and 20 concurrent users?
Vectorial data stored in PostGIS is around 20 MB only.
I need the ammount of RAM, disk space, bandwith, server cores and processor velocity.


texts: How to get detailed informations about a layer (SRID, extent, ...) in Geoserver REST API?

texts:
I'm trying go get detailed informations about a layer using the Geoserver REST API (I'm using Geoserver 2.8.1), but according to the online API documentation, I can only get basic informations such as the name and the styles.

Are these the only informations I can get using the REST API, or is there any way to get more?


texts: Thematic Graduated Bar Chart from polyline

texts:
I have a road broken into 10m sections, each section has a road condition score, I was hoping to create a bar coming of each section based on that condition and also color them based on a categories e.g 1-40 Green, 40-80 Amber, 80+ Red.

It's kind of like graduated lines but I don't want to change the thickness of the actual line, I would like for them to grow perpendicular to the line e.g

In this image each color is 10m section of road and the grey is the road. I feel this may not work well on actual road as roads curve and the height changes maybe lost but I would still like to pose the questions as it could be an inserting way of showing data along a linear feature.


texts: I want to aggregate points by count, for example, groups of 10 or 20. I use ArcMap 10.5 and ArcGIS pro, and so far i have only found tools that will aggregate by proximity and sum number of points. I need the aggregated points to be equal in weight/count.


texts: Map NDVI from MODIS image for a specific country

texts:
var dataset = ee.FeatureCollection('USDOS/LSIB/2013')
            .filterBounds(ee.Geometry.Point(-7.681, 53.43));
var image = ee.Image().float().paint(dataset, 'iso_num');

var visParams = {
 palette: ['f5ff64', 'b5ffb4', 'beeaff', 'ffc0e8', '8e8dff', 'adadad'],
 min: 0.0,
 max: 894.0,
 opacity: 0.8,
 };
 Map.addLayer(image, visParams, 'USDOS/LSIB/2013');

 // The input image has a SR-ORG:6974 (sinusoidal) projection.
 var image1 = ee.ImageCollection('MODIS/MOD13A1')
           .filter(ee.Filter.date('2016-01-01', '2016-01-31'))
           .filterBounds(dataset)
           .mean();

var ndvi = image1.select('NDVI');
var ndviVis = {
 palette: [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
'66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
'012E01', '011D01', '011301'
 ],
 };
 Map.addLayer(ndvi, ndviVis, 'NDVI');

This is my code, it works, however, the map shows the NDVI for the entire world. Is it because of the different coordination system? I tried to convert it to the same, but the Google Earth Engine API only had the projection transformation of a single image, and do not have the code for transformation of the image collection. How can I do it?

This is the link of my code https://code.earthengine.google.com/c1132b4753acf26b4ce5f93cca6ff070


texts: How do I get the extent of a rotated data frame in map units from the Layout View?

texts:
I need to get the real extent of a rotated data frame in map units from the Layout View using ArcObjects.
I've checked previous discussions on the subject:

How do I obtain the envelope of a data frame in Layout?

How can I transform layout page units to map units?

but it seems there was no final answer on how to get corner coordinates for a rotated data frame. Also, the tool provided in the above post does not recognize rotated data frames. It is a bit confusing as even if you rotate the data frame in ArcGIS, the extent, indicated in the Data Frame Properties, remains unchanged.

Is there a simple way to get the real extent directly from ArcObject, or do I need to calculate these coordinates using map scale, size of the data frame etc..? Any help would be appreciated.


texts: PostgreSQL - only allow points to be added to the table (no polygons or lines)

texts:
I want to know if there is a way to only allow certain Geometry types to be added to a table in PostgreSQL.

I have a tool that currently allows points, lines and polygons but would like to restrict it so that only points are allowed but am not sure how to do this in pgAdmin.


texts:
I'm pretty new to GIS, and since 2 months ago QGIS is a great tool to achieve very simple tasks, such open some shapefiles I find in the web and load that data on my website geoJSON formatted.

Now I'm looking forward to produce my own data, since I need the shapes of boroughs and neighbourhoods of my city and it is not available anywhere..

I've tried some things, and In my first very original idea (dream) it looked very simple... Like: open my actual city bounds layer, put and image layer behind it, resize the image to fit it, and click on that yellow pencil to fast have a beautiful hand made map... Now, after lot of hours I realize that it is not that simple...

Anyway, after some reading I could place my image on the map... I have a neigh_image_layer on which I see my full scaled beautiful JPG.

The 1st problem I faced was that my city_bounds shapefile layer was not showing anymore in the map canvas after adding the image, so.. fixing some CRS matters on image load I could place both visible, but now since my city polygon is in the right place in the middle of Brasil, my image is probably placed over an area that would wrap the whole planet.

There is any simple way to redimension the image behind the map layer? And better, am I trying the correct approach? I read tons of stuff, but maybe I don't know the correct name of the procedures I'm doing.


texts:
I'm trying to create a point shapefile from a grid index polygon shapefile I created using ArcPy at ArcGIS 10.0 . I know this is simple, but could anyone give a few pointers out there. Here is what I've done now:
but my error is:

Traceback (most recent call last):
  File "C:\Users\RLong\Documents\NewPoint.py", line 19, in <module>
    arcpy.FeaturetoPoint_management("GridIndex", "SamplesLoc", "Centroid")
AttributeError: 'module' object has no attribute 'FeaturetoPoint_management'

Any ideas?

###Simple point shapefile object for testing
###
###Import standard modules
import os, sys, math, arcpy, traceback, string, fileinput
###Setup your workspace or file you're going to use
from arcpy import env
arcpy.env.workspace = "C:/Users/RLong/GIS_Data/Conversions/GridIndex.shp"
###setup local variables
in_Features = "GridIndex.shp"
point_location = "Centroid"
out_feature_class = "SamplesLoc"

###execute point
arcpy.FeaturetoPoint_management("GridIndex", "SamplesLoc", "Centroid")


texts: How to populate field Y based on values on X field

texts: On a shapefile, how to populate fields with values based on another field. Like an if condition script. "If field X has value 1 then populate field Y with either this or that values. For example, populate values of either "Good" or "NotGood" in newly formed field based on values from field X. If field X has values 1, then populate newfield with "Good" or if the field x has value like 4 ( or any value other than 1), then populate newfield with "NotGood".


texts: How can I plot these coordinates on a map not knowing their format

texts:
Trying to plot these coordinates on a map, or convert them to Lat/Long, hoping someone can help me, way over my head. Should be somewhere in the Hess mountains Yukon Territory
y proj 890526.20462877 x proj 507373.11003899.


texts: How to create two Windows in spatstat for the same dataset?

texts:
Spatial patterns: I have a dataset of coordinates and marks (diameter) attributed to the points. I am trying to create a window (owin()) for each plot (I have 2 plots), but I don't know how to set the window for both plots:

My window:

#data: rea_spd "SpatialPointsDataFrame'
w <- owin(xrange = c(min(rea_spd$Longitude_), max(rea_spd$Longitude_)),
          yrange=c(min(rea_spd$Latitude_N), max(rea_spd$Latitude_N)))
plot(w)
#this produces a window for all the data, but I don't know how to adapt it to have a window for each plot.
#After, I calculated Kcross, which works, but I still need to have two windows, otherwise, I will violate the statistical assumptions:
kmultrea <- envelope(reav, Kcross, nsim=999, i="A", j="S",
                     global=F, savefuns=TRUE, conf.level= 0.95)
plot(kmultrea,main="Cross-K - Protected area", xlab = "Distance r [m]",
     xlim = c(0,30), ylim=c(0, 6000), legend = F)
pk <- LF.gof(kmultrea) #
pk

I tried to follow this example, but I still do not receive the points inside the windows:

m = matrix(0,20,20)
m[1+20*cbind(c$x,c$y)]=1
imask = owin(c(0,1),c(0,1),mask = t(m)==1 )
pp1 = ppp(x=c$x,y=c$y,window=imask)
plot(pp1)

Determining if trees within forest gaps are clustered using R?


texts:
I'm using mapserver and PostGis and I want to draw a geometric table in mapserver. The trouble is when I connect mapserver to PostGis and I want to draw the objects but in different colours. I mean I have this code

LAYER
NAME "poligono"
STATUS ON
CONNECTIONTYPE POSTGIS
CONNECTION "host=localhost user=*** dbname=ejemplo password=*** "
DATA " geom from prueba using srid=4258"
TYPE LINE
PROJECTION
    "init=epsg:4258"
END
CLASS
    OUTLINECOLOR 0 255 255

END

And the table has 3 columns (gid[PK], name, geom), I want to change the colour for the different names. I tried whith

CLASSITEM "name from prueba"
CLASS
    EXPRESSION "Distrito 1"
        OUTLINECOLOR 255 255 255

    EXPRESSION "Distrito 2"
        OUTLINECOLOR 0 255 255

    EXPRESSION "Distrito 1"
        OUTLINECOLOR 0 128 255

But it doesn`t work. What do I made wrong?


texts: Extract data from one shapefile based on another

texts: I am using QGIS 3.16.1.  I'm trying to find a way to have it tell me what the most common soil type is in a field boundary.  I have downloaded the soil information from NRCS's web soil survey site and have it loaded into QGIS.  In order to save time trying to eyeball what is the predominant soil type, I would like it if there was a way to have QGIS summarize either what soil type is most common or even a list of what is in each field.  I don't really know where to start, as I'm still a newbie.