texts: Why when dump a layer to memory, initial record is 1, if in original layer was 0?

texts:
I load a layer from a Shapefile (featuresA) and duplicate its contents into a memory layer (featuresB).

featuresA[0].id() is 0, but when copy his content to memory over featureB is 1,
featuresB[0].id() is 1.
Is that normal?

The following code is a  working example (standalone pyqgis script):

#!/usr/bin/env python
from qgis.core import *

layerA = QgsVectorLayer("shapes01.shp", "shapes", "ogr")

#Populate layer A
featuresA = [feat for feat in layerA.getFeatures()]

print "first id :"+str(featuresA[0].id())  # Print 0

# Duplicate layerA to layerB but into memory
layerB = QgsVectorLayer("LineString?crs=epsg:4326&index=yes", "LayerB", "memory")
attr = layerA.dataProvider().fields().toList()
layerB_data = layerB.dataProvider()
layerB_data.addAttributes(attr)
layerB.updateFields()
layerB_data.addFeatures(featuresA)
layerB.updateExtents()

# Populate layerB
featuresB = [feat for feat in layerB.getFeatures()]

print "first id :"+str(featuresB[0].id())  # ¡¡¡¡ Print 1 ¡¡¡¡¡

The code works fine and all features are copied. I would expect that the ids of the features on the two layers match, which apparently is not the case.

Why does featuresA[0].id() (from file) give me 0 and featuresB[0].id() (layer duplicated into memory) give me 1?


texts: How can I find closest point on a polygon from a point

texts:
I have a point that may be inside or outside of a polygon, I need to find the nearest point on the polygon boundary from the point.

Point(x,y) may be inside or outside of the polygon, shortest distance between polygon ABCDE and point(x,y) is marked as red line, the intersecting point is Point (m,n). I need the value of m and n.

I have ABCDE and Point(x,y)


texts: How to trigger a Python script from an ArcGIS Desktop "save" operation?

texts:
Similar to the question Run a Python Script when opening an MXD, I'm looking for a way to run a Python script that is triggered when the "Save" command/button is pressed.

I've read the Creating a Python add-in application extension help for building a Python Add-In, but I don't see an interface in the Extension class that triggers upon Save, only close.

Is there something I might be overlooking?  Otherwise, what are my alternative options other than ArcObjects, which is a road I am not willing to explore for this idea?


texts: Viewing Personal Geodatabase using software other than ArcGIS?

texts:
Is it possible to open a Personal Geodatabase file in software other than ArcGIS?

i.e are there any open source software available just for viewing the data?


texts: Decimals not working with field calculator ArcMap

texts:
I'm trying to calculate percentages of land cover in census blocks using the field calculator.

I have one field (Block_SF) that is the total square footage of the census block and another field (Imp_SF) that is the square footage of impervious surface in each block.

There is a third field (Imp_Pct) that I would like to use to calculate the percent of coverage in each block that is impervious surface.

The math is fairly simple - impervious sf divided by total square foot should yield a decimal that I can use as the percentage. However, when I use field calculator to do so, the results are always either a zero or a one, which seems to indicate that it is rounding to the nearest integer. I've tried different formats for the column - short integer and long integer - and have gone into the numeric values and made sure the decimal points were as many as four. Still the same result.

Any idea what I'm missing to get the decimal values I'm looking for?


texts: Mean monthly NDVI from Sentinel to CSV-files

texts:
I am trying to create monthly NDVI maps from Sentinel in Google Earth Engine. The script that I'm using is from things I gathered around the internet. When I try to download this to a CSV-file, I get the error: Error: Image.sample: Unable to find a crs. I know there are multiple questions regarding the monthly mean of datasets. But none try and export them to a csv-file with the coordinates and the corresponding value.
What I have so far:
    // Filter geometry and dates
    var S2_SR = ee.ImageCollection('COPERNICUS/S2_SR')
        .filterDate('2017-01-01', '2019-12-31')
        .filterBounds(roi) // filter to ROI
        .map(function(image){return image.clip(roi)}) //iterate and
                                    clips roi over whole collection
        .filterMetadata('CLOUDY_PIXEL_PERCENTAGE',"less_than", 25);

    // Add NDVI band
    var addNDVI = function(image) {
      var ndvi = image.normalizedDifference(['B8',
                                  'B4']).rename('NDVI');
      return image.addBands(ndvi);
    };

    // Apply across whole collection
    var S2_NDVI = S2_SR.map(addNDVI);

    // Select out only NDVI band
    var NDVI_only = ee.ImageCollection(S2_NDVI.select(["NDVI"],
                                                   ["NDVI"]));
   print(NDVI_only);

    //
    var months = ee.List.sequence(1, 12);
    var years = ee.List.sequence(2017, 2019);

    //calculate average values for red and nir band
    var collection = ee.ImageCollection.fromImages(
      years.map(function(y) {
        return months.map(function (m) {
          return NDVI_only
           .filter(ee.Filter.calendarRange(y, y, 'year'))
           .filter(ee.Filter.calendarRange(m, m, 'month'))
           .mean()
           .set('month', m).set('year', y);
       });
      }).flatten());

    print('done');

   //creating variables out of collection
    var listOfImages = collection.toList(collection.size());

    var s2017_1 = ee.Image(listOfImages.get(0));
    var s2017_2 = ee.Image(listOfImages.get(1));

    var sentinel2017_01 = s2017_1.sample({
      region: roi,
      geometries: true,
    });
    var sentinel2017_02 = s2017_2.sample({
      region: roi,
      geometries: true,
    });
    Export.table.toDrive({
      collection: sentinel2017_01,
      description: 'sentinel2017_01',
      folder: 'outputSentinel',
      fileFormat: 'CSV'
    });
    Export.table.toDrive({
     collection: sentinel2017_02,
     description: 'sentinel2017_02',
     folder: 'outputSentinel',
     fileFormat: 'CSV'
    });

here's the link; https://code.earthengine.google.com/b3db6273efd5a08c3089ca69dde5cae3
Anyone with any ideas?


texts: Best way to get the overlapping segment of two lines

texts:
Could anyone help me figure out the best way to get the overlapping line segment of 2 line sdo_geometry objects?

my query looks like this:

select SDO_GEOM.SDO_INTERSECTION(SHAPE, :SHAPE, 0.5) as geom
from MyTable
where MDSYS.SDO_OVERLAPS(SHAPE, :SHAPE) = 'TRUE'

the geometry objects in the table are all lines as is eth parameter being passed in, but the geometry being returned for all results is a point, not a line.


texts: How to work on histograms in ERDAS Imagine?

texts:
Is there a possibility to work on histograms in ERDAS Imagine?

For example:

Compare histograms of different images or different extents
Shift them or find patterns

I thought there must be good options as ERDAS is specialized to work on images, isn't it?

I could hardly find any histogram options in ERDAS until now, though. The only things I found was 'histogram matching' which does not leave a lot of freedom...

Are there other options to work on a histogram in ERDAS, or to observe histograms and where can I find them?


texts:
I have a file that has no extension. This file resides on an internal web server and it's called by a web map application based on django python.

It looks like :

> {"crs": {"properties": {"name": "EPSG:4326"}, "type": "name"},
> "features": [{"geometry": {"coordinates": [[127.650854900378, 1.19649
> etc...

I think it's .json, but why can't I open it? I have tried using QGIS and Arcmap, but all of them fail to open the file.

Actually, the map shows up when its opened (called) by the application, like this code below :

> var thedata = $.ajax({
>         url: thedata,
>         dataType: "**json**",
>         success: console.log("data successfully loaded."),
>         error: function(err) {
>           console.log(err)
>         }
>       });

I want to open it and export it into shapefile.


texts: How to break a large map into smaller parts (tiles), for printing, in QGIS

texts:
We have a map a very wide area. We want to break it into, for example 8, mutually exclusive parts, so we can print each one separately and, if needed, put them all together and remake the whole map.

There seems to be no way in Map Composer to automatically break the map into parts, and doing things manually would not be precise.

Is there any way to automate this?


texts: Arcgis FeatureLayer.MODE_ONDEMAND not select exactly coodinate?

texts:
I'm working with ArcGis API. I push my Feature Layer on map in MODE_ONDEMAND. But when I select feature (selectTask), the coordinates come out as [692010, 101011], with no decimal places. What I want is more accuracy, for instance: [692010.21212 , 101011.21212]. Meanwhile, in MODE_SNAPSHOT, everything is ok. Why is this?

UPDATE: Here is results MODE_ONDEMAND:

And Here is results MODE_SNAPSHOT :

Here is my code when add Feature Layer:

    var layer = vGis.Params.layers.addFeatureLayer(featureserver.url,
                { mode: FeatureLayer.MODE_ONDEMAND,
                    //mode: FeatureLayer.MODE_SNAPSHOT,
                    outFields: ["*"],
                    id: featureserver.id,
                    opacity: 1,
                }, filter);


texts: Aggregate polygons in edit session using ArcObjects and c#

texts:
I want to aggregate some polygons in edit session using ArcObjects and c# language. In other words, my purpose is to use methaheuristic algorithms for polygon aggregation. So, I want to aggregate some polygons, see the numerical results and if it is not good, don't save it, and aggregate some others. And, in order for that, I want to implement it using edit session.For that purpose, I have implemented "Merge" which is in the editor toolbar. The code that I have used for merging two features is as follows:

Private void MergeFeatures (IFeature mainFeature, IFeature subFeature)
{
ITopologicalOperator union=mainFeature.Shape as ITopologicalOperator;
mainFeature.Shape=union.Union (subFeature.Shape);
mainFeature.Store();
}

But my problem is that it does not bring the polygons together. In other words, it just changes the attribute table of the feature class. The screenshot of part of my data is as follows:


texts: Python collect all points that meet a criteria for every point in the shapefile

texts:
I have a shapefile of points that has a field of yield data from harvesting. I am trying to 'smooth' the data by taking each point and collecting all points within 50ft and then update the shapefile with the average of the 50ft. circle. The end goal is that I want a shapefile with every point having a field that contains the average yield of all the other points within 50ft.

Unfortunately python isn't my main programming language and I'm just learning about gis.

The code below is working, but I'm not satisfied with it. What is a performant way to find all points within a given radius of another point, and then do that for every point in the file?

The outliers_r is a field on the shapefile where I cleaned up some of the outliers in the data.

The yield_smoo is the field I want updated with the average yield.

from osgeo import ogr
import fiona
import pyproj as proj
from shapely.geometry import shape
from shapely import geometry

#open the file
source = ogr.Open("2017_wheat_yield_data.shp", update=True)
layer = source.GetLayer()

fc = fiona.open("2017_wheat_yield_data.shp")
fc_2 = fiona.open("2017_wheat_yield_data.shp")

crs_wgs = proj.Proj(init='epsg:4326')
crs_bng = proj.Proj(init='epsg:3420')

for feat in fc:
    coords = feat["geometry"]["coordinates"]
    x, y = coords
    x1, y1 = proj.transform(crs_wgs, crs_bng, x, y)
    point_1 = geometry.Point(x1, y1)

    yield_data = []

    for feat_2 in fc_2:
        coords_2 = feat_2["geometry"]["coordinates"]
        x_2, y_2 = coords_2
        x2, y2 = proj.transform(crs_wgs, crs_bng, x_2, y_2)
        point_2 = geometry.Point(x2, y2)

        distance = 50
        # circle_buffer = point_1.buffer(distance)
        # if circle_buffer.contains(point_2):
        #     print('circle buffer contains point 2')

        if point_1.distance(point_2) < distance:
            data = feat_2["properties"]["outliers_r"]
            yield_data.append(data)

    # sum the yields
    average_yield = sum(yield_data) / float(len(yield_data))
    id = int(feat["id"])
    fture = layer.GetFeature(id)
    fture.SetField("yield_smoo", average_yield)
    layer.SetFeature(fture)


texts: QuickMapServices not showing basemap for my location - QGIS

texts: I've plotted several locations based in Auckland, New Zealand, using Grid NZGD2000:ESPG:2193. When I load QuickMapServices and set to the same grid, The basemap produced is situated off Omen.


texts: How to find the attributes of a shapefile from the command line?

texts: I'm following Mike Bostock's tutorial "Let's Make A Map." He explains how to filter a shapefile using ogr2ogr based on the attributes of the features. I can open a shapefile in QGIS and look at the attribute table to find out what those features are. Is there a way to find that information from the command line?


texts: Creating shapefile from lat/long values using ArcPy?

texts:
How can I create a shapefile using Python in ArcGIS 10?

I have lat & long.

From this I need Python code which will create a shapefile in ArcGIS Desktop 10.


texts: Publishing ArcGISTiledMapServiceLayer on ArcGIS for Server?

texts:
Please suggest some URL or steps to publish ArcGISTiledMapServiceLayer service on ArcGIS 10.1 server so I can add in my application.

What are the best practice to do that to increase performance?


texts: Proj4 string missing datum information for SAGA outputs loaded into R

texts:
I've produced a set of topographic derivatives using both SAGA and ArcGIS. They are all in the same CRS.

When I load them into R, however, the proj4 strings are different:

>arcRaster
class      : RasterLayer
dimensions : 16002, 21200, 339242400  (nrow, ncol, ncell)
resolution : 10, 10  (x, y)
extent     : 480200, 692200, 4896280, 5056300  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

> sagaRaster
class      : RasterLayer
dimensions : 16002, 21200, 339242400  (nrow, ncol, ncell)
resolution : 10, 10  (x, y)
extent     : 480200, 692200, 4896280, 5056300  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Note the raster produced in saga is missing +datum=NAD83. They align properly and they stack without problems. Why is the datum information missing?


texts: What does "Remote Styles Cache" option of GeoServer do?

texts:
There is no information about the parameter in the GeoServer documentation.

Can you explain what happens when the parameter is activated and will it speed up the work with a large number of layers?


texts: Transforming multiple images from projection EPSG:3857 to EPSG:21781

texts:
I have several images in EPSG:3857 which are originally from the Google maps API. I have a map layer which is based on OpenLayers v3.6 and it is in EPSG:21781. I would like to project the images from EPSG:3857 on top of EPSG:21781 layer. From the Openlayers API V3.11 there is an option to reproject at the client side on the fly with the following example code.

layers['bng'] = new ol.layer.Tile({
        source: new ol.source.XYZ({
          projection: 'EPSG:3857',
          url: 'https://tileserver.maptiler.com/miniscale/{z}/{x}/{y}.png',
          crossOrigin: '',
          maxZoom: 6
        })
      });

Since I am stuck with V3.6, I would like to convert the images at the server side and project them on top of EPSG:21781. Could someone help me with this.


texts: In ArcGIS, I have a grid 5x5km (I have it both in vector and raster formats) I want to calculate, based on this grid, the zonal statistic (MEAN) of a raster 30x30m pixel size. The raster overlaps perfectly the reference grids, the output (which I set in the environments/Raster_analysis to be 5X5km) is somehow shifted of about half a grid (which does not allow me to transfer the MEAN data from the raster to the feature grid, throughout RasterToPoint and then join by location).


texts:
I have Microsoft MapPoint 2006 installed as a part of a larger system for home/social care planning. The maps and the address information points are quite incomplete and I would like to roll out my own layers that are way better.

Can anyone explain how to do this or perhaps point to where I can find the information on the subject?


texts: Creating markers along a border in Google Maps

texts:
Is it possible to create markers along a country's border in Google Maps?

All the posts I've been able to find only give advice on how to create markers along the roads.


texts: Extracting trees/vegetation as separate point cloud using ArcGIS Pro

texts:
I have ArcGIS Pro 2.4 with 3D Analyst and Spatial Analyst Extensions. I would like to extract all points from a point cloud (.las) file that match a certain classification to a new, separate point cloud.
what I have tried:
I have ran the "Classify LAS By Height" tool and have assigned a class to points matching a certain height above ground. When I filter the point cloud by this new class (change symbology) I see the results I expect, which is great, but I want to extract/save-out these newly "classified" points to a new point cloud (.las) file.
How do I do this?


texts:
I am using openlayers3 and i am initializing my map like so:

var view = new ol.View({
    center: ol.proj.fromLonLat([5.8713, 45.6452]),
    zoom: 12,
});

var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM()
        })
    ],
    target: 'map',
    controls: ol.control.defaults({
        attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
            collapsible: false
        })
    }),
    view: view
});

Now how do i get the extent array , i tried the following code:

var extent = map.getView().getExtent();
console.log(extent);

But that does't seem to be the right code, i get an erros saying

  Uncaught TypeError: map.getView(...).getExtent is not a function

can anybody explain , how do i go about getting the correct extent ?


texts: Free satellite data availability for moon

texts: What are the satellite raster and vector data sets available for mapping the moon that can be freely downloaded? I need to know the availability and type of the data sets from which mineral mapping, feature mapping and other analysis of lunar surface can be performed.


texts: Transforming and Projecting Shapefiles in R & ArcGIS for Desktop?

texts:
Perhaps this is an issue with one of my original data sources, but I've searched extensively to address this problem both through ArcGIS and R.

I've got 2 separate sources of a county layer SourceA and SourceB.

I'm trying to project the SourceA layer so that it aligns with the SourceB layer.

The original SourceA layer is in NAD83, UTM zone 10 projection:

proj4string :[+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs+ellps=GRS80+towgs84=0,0,0]

The SourceB layer is in NAD83, Albers, California.

proj4string : [+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs+ellps=GRS80 +towgs84=0,0,0]

I've tried to transform the county data using both R and ArcGIS, but am having problems.

When I try to project it in ArcMap, the program puts in a string of characters (see circled image below: WGS84 to NAD 83 + WGS84 to NAD83) in the geographic transformation section, which it shouldn't because they're both supposedly in NAD83... (right?), and it won't let me take that string off the chart with that little X button.  If I run it anyways, as setup in the screenshot below, it throws an error "invalid extent for output coordinate system".

When I use use spTransform in R to project the SourceA data to Albers California, it runs the transformation and results in a skewed county shape in comparison with the original SourceA shapefile (which is in UTM). I think this is to be expected (see screenshot below).  But... that transformed layer won't align or even show up on the same map with the SourceB data that is also projected in NAD83, Albers California .

The SourceA layer after spTransform to Albers California reads the following proj4string:

proj4string :[+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0+y_0=-4000000 +datum=NAD83 +units=m +no_defs+ellps=GRS80 +towgs84=0,0,0]

Here's the R code using the spTransform command to transform the original SourceA data into NAD83 Albers, California resulting in the image above. Sourced from spatialreference.org:

county <- spTransform(county,CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

Three zipped shapefiles illustrating the problem can be downloaded off this dropbox link.  1) SourceA original projection, NAD83 UTM zone 10, 2) SourceA transformed to NAD83 Albers California through R's  spTransform, and 3) SourceB in NAD83 Albers California.

Is there something wrong with the original SourceA or SourceB projection information? Any suggestions?


texts: Changing raster calculation values - shorten this script?

texts:
i'm doing something quite straightforward regarding raster algebra but am struggling either to find the correct function or make a function work correctly;

Basically I have 2 rasters, representing consecutive years, that represent a classification (including NA values);

r <- raster(ncol=10,nrow=10)
r[] <- sample(c(1,2,4,8),size=100,replace=T)
r[runif(10*10) >= 0.50] <- NA

r1 <- raster(ncol=10,nrow=10)
r1[] <- sample(c(1,2,4,8),size=100,replace=T)
r1[runif(10*10) >= 0.50] <- NA

to get my change from one year (r) to the next (r1) I simply subtract r1 from r;

r2 = r - r1
freq(r2,digits=2)

Each new integer codes to a specific change, 0 simply means no change. I don't mind that cells become NA if NA only exists in one raster layer, this is fine. What i do want to do is to examine the 0 values more closely so if the original value in r is 1 and the recent value in r1 is also 1, i want to know this in the resulting r2 - ie quanitfy and analyse the 'no change' cells. Same for r=2 & r1=2, r=4 & r1=4 etc  - not just all as 0s but as a collection of new, unused codes to represent what cells stayed the same (and how they stayed the same) from one year to the next.

Basically, i extracted the position of the cells from r2 that equal 0 as a spatial points data frame;

pts <-rasterToPoints(r2,fun=function(x){x==0},spatial=T)
plot(pts)

then i replaced the cells in r2 that equal 0 with values from r (or r1, doesnt matter as they are the same value) using the locations derived above;

r2[r2==0] <- (extract(r,pts))/10
plot(r2)
freq(r2,digits=2)

dividing the replacement values by 10 ensures they do not fall into the same 'bin' as any other mapped code change.

I am sure there is a quicker way? i thought i'd be able to create a raster stack from r and r1, creating an r2, then use some sort of 'where' or 'ifelse' function to perform the calculations above, however everything i do either results in errors or neglect of NA values to make the functions work. Or S4/integer errors.

I'm still not sure this will work on a big sample yet, and my final datasets are very massive.


texts: Updating geomesa record using Geomesa-Accumulo and Spark with Java API?

texts:
I have found here how to save data using spark and geomesa accumulo. But I can't find anywhere how to update existing records.

I know that Accumulo is under the hood. Accumulo provides data versioning (see this link). Do I have to use versioning of data? I don't know where to find any geomesa documentation concerning records' update.


texts: Trying to install gdal 3.0.4 on Red Hat 8

texts:
I'm running Red Hat 8, and I'm trying to install gdal 3.0.4 on it.

The commands I've excuted thus far:

wget http://download.osgeo.org/gdal/3.0.4/gdal-3.0.4.tar.gz
tar zxvf gdal-3.0.4.tar.gz
cd gdal-3.0.4
./configure
sudo make install

After I run sudo make install, I get the following error messages:

path/to/gdal-3.0.4/.libs/libgdal.so: undefined reference to `jpeg_CreateDecompress_12`
path/to/gdal-3.0.4/.libs/libgdal.so: undefined reference to `jpeg_write_tables_12`
path/to/gdal-3.0.4/.libs/libgdal.so: undefined reference to `jpeg_start_decompress_12`

...

It continues on like this for some amount of lines all starting with jpeg and ending with 12 and different words in between. It finally ends with:

collect2: error: ld returned 1 exit status
make[1]: *** [GNUmakefile:82: gdalinfo] Error 1
make: *** [GNUmakefile:112: apps-target] Error 2

I'm not sure what to do or what this error even means. If anyone has any ideas, or if there is anymore information I can share, let me know.


texts: Writing an if then statement for if a parameter is selected (not Null) in Model Builder

texts:
I have a tool I'm trying to create in model builder. I have run into a problem as I have basically no experience with python. I'm trying to basically use the calculate value tool and an if then statement to determine if a model parameter (type String) has a value inputted to it by the user, if it does then go one direction in model builder and if not go another. Something to the effect of the picture below:

The Second picture is of the model with the properties of the string and the tool (in the lower right)


texts: How do I get the upstream distance to all points in a DEM from a single downstream point?

texts:
I want to get the upstream distance to all points in a DEM from a single downstream point. The goal is to create a "downstream distance to a point grid". This then would allow me to create "distance" contours upstream from a point in the watershed.

I tried using the Raster Calculator with this equation:

Outras = Con("Fdr" == 1,16,Con("Fdr" == 2,32,Con("Fdr" == 4,64,Con("Fdr" == 8,128,Con("Fdr" == 16,1,Con("Fdr" == 16,1,Con("Fdr" == 32,2,Con("Fdr" == 64,4,Con("Fdr" == 128,8)))))))))

This "reverses" the flow direction in the Fdr grid, but doesn't really work.

Any suggestions?


texts: Labelling features by multiple attributes in QGIS

texts:
I am a fairly competent QGIS (3.14) user but really struggling with a particular labelling issue. I want to create attribute-based symbols/labels based on multiple attributes for a single point, but with different combinations based on attribute values.
My data is structured as follows: each geographic point is a row, and then I have a series of attributes (hiking, swimming, fishing) as boolean yes/no values. What I want is to create labels or symbols for each point, and custom labels/symbology based on the boolean values (so, if place A has fishing, show a fish; for place B with fishing and swimming, show a fish AND a swimmer icon).
This would like something like this:

but with variation in the "Hiker/Biker Camp" image (not my map - an example to share).
I've tried using rule-based symbology with manual offsets, but that leaves an ugly gap if a point has attribute A and C but not B.
Is there a simple way to do this in QGIS 3.14 or with plug-ins?


texts: How to convert from .shp to .csv using python and either gdal_translate of ogr2ogr

texts:
I am trying to convert .shp files to csv using either gdal_translate or ogr2ogr and python but its not work. if i run the following from qgis its not giving an error but at the same time its not producing the output files. Please help. Below is what  i was running:

import os
command = "ogr2ogr -f 'CSV'" + "C:\Users\zanamwe\Desktop\Desktop\Programming For GIS\shapefiles\lakes.csv"+ "C:\Users\zanamwe\Desktop\Desktop\Programming For GIS\shapefiles\lakes.shp"
os.system(command)


texts: Pyqgis 3: Trying to Update Attribute with Value but Getting Index out of Bounds

texts:
I'm able to successfully add attributes with an expression to my vector layer using:

#Add fields to get evaulated based on an expression
vectorLayer.addExpressionField("'R'",
QgsField("rtype",QVariant.String, 'string', 11))

vectorLayer.addExpressionField("y(point_on_surface($geometry))",
QgsField("lat", QVariant.Double, 'double', 20, 7))

vectorLayer.addExpressionField("x(point_on_surface($geometry))",
QgsField("long", QVariant.Double, 'double', 20, 7))

vectorLayer.addExpressionField("\'Comment explanation; type: {}\'".format(type), QgsField("comment",QVariant.String))

vectorLayer.updateFields()

Then based on a select by location I need to update the comment value, so I get the indexes using:

request = QgsFeatureRequest().setFilterFid(1)
features = vectorLayer.getFeatures(request)

#Get indexes for attributes we are keeping
for feature in features:

     rTypeIdx = feature.fields().indexFromName("rtype")

     latIdx = feature.fields().indexFromName("lat")

     longIdx = feature.fields().indexFromName("long")

     commentIdx = feature.fields().indexFromName("comment")

To update the comment accordingly I use:

selectedFeatures = vectorLayer.selectedFeatures()

for feature in selectedFeatures:

     fid = feature.id()

     #Get comment index and value for updating attribute values
     commentIdx = feature.fields().indexFromName("comment")

     commentValue = 'Comment value'

     attributes = {commentIdx: commentValue}

     feedback.pushInfo("Updating feature {} with attributes{}".format(fid,attributes))
     vectorLayer.dataProvider().changeAttributeValues({fid: attributes})

After I run my script everything works, except I get the warning that my attribute index for the comment is 10 out of bounds [0;7].  I've verified this index using several different methods including feature.fields().indexFromName("comment").  They all return the index value of 10 which is the last attribute.  If I hardcode the value of 7 for my index it crashes QGIS.  What am I doing wrong?  How can I get the correct index value to update this field without getting an index error? I'm using QGIS version 3.x.

EDIT:  I'm now trying to update my feature using

vectorLayer.startEditing()
featureUpdated = vectorLayer.changeAttributeValue(fid, commentIdx, commentValue)
vectorLayer.commitChanges()

This is also failing by returning False.


texts: Mosaic Aerial Photography Tiles to a Seamless Raster

texts:
I've got about 6 different Aerial Photography datasets I am managing at the moment. All of them are made up of tiles covering a large area. Some of these datasets contain up to 2000 tiles. All the tiles are slightly overlapping so there re no gaps in the data. Tiles work fine for me, but other users are requesting Mosaic's of the datasets.

I've been trying over the last few days to Mosaic these tiles into a seamless raster. I've been using the ArcGIS Mosaic to New Raster Tool. It worked great on the first (smallest) dataset, but since then ArcGIS crashes after running for about 8 hours. I've tried the tool on a more power machine and with less tiles, but ArcGIS still keeps crashing.

I'm not sure what else to do here. I need to mosaic these tiles into one raster, but despite any changes to the options I keep getting the same result. I'm wondering if the Moasic to New Raster is the correct tool for this job, or should I be using something else?

I'm looking for a either a ArcGIS or python programming solution to this issue.

The options I am using for the Mosaic to New Raster Tool:

Input: 500 to 2000 input tiles that are JP2
Output: Save the new raster as a JP2
Coordinate System: NAD83 17N to NAD83 17N
Pixel Type: 8 Bit Unsigned to 8 Bit Unsigned
Cellsize: 0.2 to 0.2 (varies, but the the cellsize stays the same from the input to the output)
Band: 3 (4 in two dataets)
Mosaic Method: I've tried pretty much every option here, the preferred is Blend.
Colourmap: First (I haven't changed this at all yet)


texts: How to add rings that intersect with a polygons edge in QGIS

texts:
I need to amend the area of an existing polygon so that I can draw another polygon (of another typology in the amended area). Normally, I would manually move the nodes but this will be very time consuming due to the curving edges of the polygons.

Is there a way I can use the add rings button to do this?

(so far it is not letting me amend an area that intersects the edge of the polygon)

Or if not, is there an alternative option?

The area I need to remove lies above the footpath in the screen grab below


texts: Canadian Shapefile from Statscan is wrong?

texts:
It's great that Stats Can produced this for free!  However, it's useless if it's wrong.  I'm using the shapefile for territory planning so the FSA (postal code) data needs to be 100% accurate.  But look at this image!  Why would there be the same FSA all over the place?

Where can I find a 100% accurate shapefile for Canada?  I'm fine with paying if I need to...

Edit: Here's another one from my old home postal code T5E.  Why the heck is the Statscan shapefile showing T5E above Edmonton?  If you look at Google Maps, it's got it right...


texts: Is there any way to include "universal" style/symbology in GeoPackage?

texts:
I am just playing with the GeoPackage format in QGIS. In QGIS it works nice with the style saved from QGIS but if I load it for example in SAGA the style is not applied.

Is there any way to include at least some basic style (for example point colors for value ranges or classes) so it could be used by different applications?
Or is the included style always "application-specific"?


texts: PostGIS ST_Intersects vs ArcGIS Select by Location

texts:
I have two polygon tables, 'cities' with ~5600 rows and 'lte_coverage' with ~469314 rows. I need to know wich cities have LTE Coverage.

Using ArcMap, I get this answer in about a minute, using 'Select By Location', but using PostGIS it takes really long, indeed after 90 minutes I cancelled the query. I am wondering what can I do to get a better performance on PostGIS.

Table cities:

                                                        Table "public.cities"
   Column   |          Type           |                           Modifiers                            | Storage  | Stats target | Description
------------+-------------------------+----------------------------------------------------------------+----------+--------------+-------------
 id         | integer                 | not null default nextval('cities_gid_seq'::regclass)           | plain    |              |
 geom       | geometry(MultiPolygonM) |                                                                | main     |              |

Indexes:
    "cities_pkey" PRIMARY KEY, btree (id)
    "cities_geom_idx" gist (geom)

Table lte_coverage:

                                                        Table "public.lte_coverage"
   Column   |          Type          |                              Modifiers                              | Storage | Stats target | Description
------------+------------------------+---------------------------------------------------------------------+---------+--------------+-------------
 gid        | integer                | not null default nextval('lte_coverage_gid_seq'::regclass)          | plain   |              |
 geom       | geometry(Polygon,4326) |                                                                     | main    |              |

Indexes:
    "lte_coverage_pkey" PRIMARY KEY, btree (id)
    "lte_coverage_geom_idx" gist (geom)

I am using this query

SELECT distinct c.id FROM cities c, lte_coverage l WHERE ST_Intersects(c.geom,l.geom)

And I have this explain output:

HashAggregate  (cost=304837.32..304892.99 rows=5567 width=8)
  Group Key: c.id
  ->  Nested Loop  (cost=0.29..302660.09 rows=870890 width=8)
        ->  Seq Scan on cities c  (cost=0.00..4715.67 rows=5567 width=40)
        ->  Index Scan using lte_coverage_geom_idx on lte_coverage l  (cost=0.29..53.36 rows=16 width=1132)
              Index Cond: (c.geom && geom)
              Filter: _st_intersects(c.geom, geom)

I am using

gisdb01=# select version();
                                                 version
---------------------------------------------------------------------------------------------------------
 PostgreSQL 9.5.5 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-4), 64-bit
(1 row)

gisdb01=# select postgis_full_version();
                                                                            postgis_full_version

------------------------------------------------------------------------------------------------------------------------------------------------------------------------
----
 POSTGIS="2.2.4 r15258" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" TOPOLOGY RAS
TER
(1 row)

So, is there anything I can do  to PostGIS perform better?


texts: Rasterizing all intersected cells when rasterizing line layer in QGIS

texts:
I have a line layer that I want to rasterize. I run the Rasterize (vector to raster) tool in QGIS Processing to do so. However, the output is lacking in that it doesn't capture each cell that the line layer touches. See image below.

So, how can I rasterize a line layer in QGIS and make sure all intersected cells are given a value?


texts: Using Python to select records by date field

texts:
I am trying to write a python script to select records created within the past 7 hours (date field - database time) using a where-clause within arcpy. The table is ArcSDE database (Microsoft SQL Server).

I've written a script that calculates the time 7 hours prior to the current time. Then I want to use the "report_time" variable to select the relevant records in a table view but I'm getting an invalid expression error. I've tried to reformat the SQL statement in every way I could think of and still get the invalid expression error (I've kept them in my script and commented them out for reference).

-- Get start and end times for report

start_time= datetime.timedelta(hours = 7)
end_time = datetime.datetime.now()

report_time = end_time-start_time #this is the time that gets used to filter records

-- find all records that are later than or = to report_time

SQL = "created_date >= " + report_time.strftime('%m/%d/%Y')
print SQL

arcpy.SelectLayerByAttribute_management(ViewTable,"NEW_SELECTION", SQL)

-- SQL = "created_date >= " + str(report_time.strftime('%m/%d/%Y'))

-- SQL = '"created_date"<='+ report_time.strftime('%Y/%m/%d')

-- SQL = '"created_date"<='+ report_time.strftime('%m/%d/%Y')

-- SQL = "'created_date'<= "+ report_time.strftime('%m/%d/%Y')

-- SQL = '"created_date"<='+ report_time.strftime('%m/%d/%Y %I:%M:%S %p') - ExecuteError: ERROR 000358: Invalid expression

-- SQL = "'created_date'<= "+ report_time.strftime('%m/%d/%Y %I:%M:%S %p') - ExecuteError: ERROR 000358: Invalid expression

-- SQL = "'created_date'<= "+ str(report_time.strftime('%m/%d/%Y %I:%M:%S %p'))

-- SQL = '"created_date"<='+ str(report_time.strftime('%m/%d/%Y %I:%M:%S %p'))

--SQL = '"created_date"<= report_time.strftime'('%m/%d/%Y %I:%M:%S %p') - TypeError: 'str' object is not callable

--SQL = '"created_date" <= report_time' #this returns an expression error

--SQL = "'created_date' <= report_time" #also tried this - expression error

--SQL = 'created_date'<= report_time # returns error: TypeError: can't compare datetime.datetime to str


texts: Tracing Upstream Accumulation from a polygon start

texts:
I am new to ESRI products (I typically use QGIS), but am using it in order to analyze upstream/downstream flowlines with the National Hydrography Dataset (specifically, HUC-4 mapping units for NHDPLUS HD data). In particular, I am hoping to select all upstream and then downstream flowlines from a polygon layer. An example might be - what are the upstream flowlines that flow into this state park? It seems like the only way to do this is the trace upstream flow accumulation by point (by dropping a flag). This will be prohibitively cumbersome to try and drop a flag at all of the points the polygon boundary intersects with NHD flowlines.
Is there a way to set the starting point from the polygon object itself? I have state parks in one layer, and NHD flowlines in another. My ultimate goal is to be able to identify upstream/downstream monitoring stations (which is a point layer of water quality monitoring stations).


texts: Cloud masking with the Descartes Labs Platform Scenes API

texts: I have created a scene collection in the Descartes Labs Platform using the Scenes search function (dl.scenes.search) for my area of interest and sensor of choice (Sentinel-2). I would like to iterate through the scenes in the collection and mask for clouds using the cloud-mask band before creating a composite with only clear pixels. How would you approach this task?


texts: Creating polygons from one to many relationship of points using QGIS or PostGIS?

texts:
I have a set of polygons that are in a one to many relationship with points.

The polygons are administrative boundaries and the points have information on companies and percentage of ownership .

I need to get a summary of the company data for each polygon.

The way I know how to do it is to do a relate. I did a spatial join on the points to attach the polygon information and then did a relate on the polygon table to highlight points that fall within it. From there I can do a summarize to get the values I need (sums/averages of company data within each polygon).

But that process is only good for one at a time processing. I have several hundred polygons I want the information for.

Any ideas?


texts: Background label shows above all annotations and layers

texts:
I'm using ArcMap 10, and have a very dense map of Detroit with many point features.  To lay out the labels properly, I converted the labels to annotation (stored in the map, not in a geodatabase), and spent hours and hours arranging them.

Now I'm trying to place a background label for Detroit.  I'm using the Maplex engine, and in the label properties, I've selected "Background label (placed first)".  But the label is still on top of everything else -- all other layers and the annotations.  I also tried checking "Never remove (allow overlap)".  How can I put it in the background?  Also, are annotations always drawn above all labels?  It seems like that's true.


texts:
I have a remote machine where my GRASS binaries (7.8) are available and I want to apply this solution https://gis.stackexchange.com/a/255899/36712 to my data set of thousands of small images. In most PyGRASS tutorials I find the API works more like an remote to the GUI, in other words you apply functions to layers present in the GrassGUI for example gscript.parse_command('r.univar', map='elevation', flags='g') where I apply the function r.univar to the layer elevation
I allready found the grass-session libary allowing me to create a session without "starting" GRASS. Most of the tutorials I find are also either old referring to Python 2.7 or don't realy cover how to load a dataset or what would even better how I could pass a variable containing the file object to GrassGIS. How I imaging the code should look like is this:

from grass_session import Session

with Session(gisdb="/tmp", location="location", create_opts="EPSG:4326"):
    for png in list_of_png_files:
        geo_tiff = transform_png_to_geotiff_with_gdal(png)

        #either i need to do this:
        Loaded_geotiff = load_geotiff_with_grass(geo_tiff)

        #or I can work directly with the geo_tiff Variable
        output = processing.runalg('grass:r.mapcalculator',geo_tiff)
        thined = processing.runalg('grass7:r.thin',output)
        lines = processing.runalg('grass7:r.to.vect',thined)

        save_lines_to_shp(lines)

How can I achieve this? The PyGRASS API seems so counterintuitive and it's hard to understand how this is supposed to work.


texts: Cache lookup failed for function 1 while activating PostGIS

texts:
I'm using PostgreSQL 10.5 on a remote server (Debian). I created a database inside then I attempted to enable PostGIS but it wasn't successful.

I had this error : cache lookup failed for function 1


texts: Batch checking of 1000 records online without batch option. Need automation

texts:
Do you know of any method for looping 1000 records over this website?

https://apidocs.os.uk/docs/os-places-match

I would like to check some addresses but OS doesn't provide batch checking.


texts: Joining multiple fields to table and export new table as join field name using ArcGIS Geoprocessing?

texts:
I am using join field tool to join table (info tables/.dbf) to one  shapefile.

I am trying to use ModelBuilder, as there are around 100 tables (each with 5 columns X 50-150 rows), each has to be joined with one shapefile (3 column X ~2 million rows). After the join new shapefile has to be saved as respective field table’s name.


texts: I've been using arc.js and I was trying to see if there was any way to color the arcs, except for the normal blue shown(ex). I tried looking around with the arc options, but I couldn't really find anything too useful. Is there anyway I can do this with arc.js?


texts: Is it possible to convert QGIS QLM to SLD (w/ New Symbology)?

texts:
I have a QGIS project where I have created some nice styles using New Symbology. I would like to convert this to either SLD for use in GeoServer or CartoCSS for use in TileMill. I have tried saving the style as SLD directly from QGIS and a number of plugins, but I always get "New Symbology not supported" or "<!--FeatureRendererV2 VectorFieldRenderer not implemented yet-->".

I have also tried a number of plugins to convert the project/layer to tilemill or a mapserver project, but no dice.

Any help with this would be appreciated.


texts: Importing OpenStreetMap (OSM) to ArcMap?

texts:
I am having some problems with the import of a locally saved osm file into arcmap 10.1. I have installed the Arcgis Editor for osm from here: http://www.arcgis.com/home/item.html?id=6a2a3c3cece749558393d4e80241ef51.
I am using windows 10 64 bit. The arcgis 10.1 seems to work fine.
I am trying with Download Extract Symbolize tool and this is the message it gives me.

There we go this is what happens when i try to load a locally saved OSM file with JOSM editor.


texts:
I'm loading OpenStreetMap data into a PostGIS database with Imposm3.

Regarding the documentation Documentation, fields with NULL are stored as empty field.

To query attributes who are not empty I try:

SELECT * FROM table.polygon WHERE 'attribute' is not ''

Unfortunately, this is not working.

How do I exclude empty string fields as query in a PostGIS database?


texts: FlowMapper plugin Error in QGIS 2.12.0 Lyon

texts:
I have a problem using FlowMapper plugin. But I don't now why. This is the message:

Traceback (most recent call last):
File ".../python/plugins\FlowMapper\flowmapper.py", line 2656, in run
flowpyv07.shapefilemaker(FlowType,CreateShpNodes,IncludeNodeNames,str(SaveDirectory), str(SaveShpName),str(SaveShpNameNodes),str(InputMatrixName),str(InputNodesName),str(InputNodeNamesName),str(combotext))

File ".../python/plugins\FlowMapper\flowpyv07.py", line 51, in shapefilemaker
onevalue = float(separatestrings[columns])
IndexError: list index out of range

Python-Version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS-Version: 2.12.0-Lyon Lyon, cd9d645


texts: Rasterize/Filter with Burn and Weight variables in QGIS

texts:
Using QGIS 3.10.1, I am attempting to rasterize some demographic polygons, and then apply a filter (Gaussian) to create a smooth/graduated surface. This works fine for some of my values (e.g. total population and housing units), but I am running into problems rasterizing averages (e.g. age and income). Average age, median income, and total population are all attributes of each demographic polygon.

When I rasterize the average age, for example, the values in the raster output are not weighted to the number of people, and this error gets magnified when I apply a filter. Imagine a feature that has 1000 people and an average age of 60 and the surrounding polygons have 10 people each and an average age of 20. When you use rasterize it creates a cell that is bigger than and contains all these features, but the rasterize algorithm only asks for the burn value and as a result, it averages the average age of all the features without regard to the number of people. As a result, the raster cell for that area erroneously shows a value of 22. If it was correctly weighted for the number of people it's value would be 59.

Because the area that I am looking at is a large area and includes urban polygons that are very small and rural ones that are very big it is inevitable that some cells span multiple polygons, and that there will be some overlap.

Are there any tools that can rasterize and apply a filter using 2 variables (a burn variable and a weight/quantity variable or can anyone recommend a workflow to address this issue?

One possible solution I've been exploring is to create duplicate geometries based on the total number of people. For example, a polygon with an average age of 40 with 10 people would get 1 polygon and a neighboring polygon with an average age 35 over 100 people would get 10 identical polygons so that when they go into the rasterize algorithm they are weighted appropriately. That could solve the rasterize issue, but I still would have the weighting problem when I apply the Gaussian filter.

Another way I am thinking about solving this issue is to rasterize/filter using aggregate data, and then use the raster calculator as the final step to divide the aggregate data by the total numbers of population or households.


texts: color polygons based on 2nd layer lat and long

texts:
I have nyc PLUTO map as layer one (parcel map showing polygons of building outlines)
I have a second layer with Lat and Long identifying categories of buildings.

I would like my layer one polygons to be colored if they have the second layer intersecting with them


texts: Stacked Label using If/Else in a Label Expression

texts:
I have a large set of data with lab results that I'm looking to label in ArcMap. I'm trying to build a label expression that needs to produce a label that essentially needs to look like this:

B=5
T=10
E=200

So if B [BEN] is >= to 5, it will produce text that is bolded and red. If B is < 5, it will come out green, and if it is anything else it will come out B=ND. This is the code I have so far to do that, and it works:

def FindLabel ([BEN]):

 if float([BEN]) >= 5:
  return "<CLR red='255'><FNT size = '8'><BOL>" + "B=" + [BEN] + "</BOL></FNT></CLR>"

 elif float([BEN]) < 5:
  return"<CLR green='255'><FNT size = '8'>" + "B=" +  [BEN] + "</FNT></CLR>"

 else:
  return  "<CLR green='255'><FNT size = '8'>" + "B=ND" + "</FNT></CLR>"

The problem is that I also need to do this for compounds T, E, and so on...If I didn't want to use the if/else it would be simple, I could use something like:

return "B=" + [BEN] + "\n" + "E=" + [ETH] + "\n" + "T" + [TOL] + "\n" + "I=" +  [ISO]

I basically need to have stacked labels (as shown above) while using the If/Then statements, and I'm really blanking on how to do this.


texts: SAGA Gaussian filter lower value than 1?

texts:
I'm smoothing a 2 meter DEM using the Gaussian filter in SAGA through QGIS/Processing and the smallest amount for the smoothing of the Standard Deviation value of 1 which for me is too much, it all gets way too blurred. Is there another way of smoothing my DEM with either another approach (I use Mac) or tweaking SAGA so that the smoothing will be less?

I guess one way would be to resample my data to 1 m DEM and then run the filter again?

DEM

DEM with gaussian SD 1


texts: Symbolizing lines based average value of point (start, middle, end) in ArcGIS Desktop?

texts:
There are three points in each line of road network, each point has weight and I want to give the average weight of these three points to corresponding line and symbolized all the lines accordingly. Any idea how to do it?! To remind you both points and lines are two different shape files.
Point values are centrality measures!


texts: Write SearchCursor records into a table in reportlab

texts:
I started to use ArcGIS Pro 2.0 recently. There is no reporting yet and I am using reportlab to create my pdfs (.rlf does not exist in ArcGIS pro). I also use platypus to create tables. The forums are very active so this part was okay.

I am working with parcels. I have to export the attribute tables records, after making a selection of the new parcels created or updated (landsurveyors need to review them). I am using arcpy.searchCursor to open and read my table. I can print the results after sorting them. However I would like to write the results in a table format Inside my pdf template.

I found this: Using ArcPy and Python to export Table View as PDF format?, which is close to what I want. However this is for ArcMap and does not seem to work with ArcGIS pro.

Any ideas?


texts:
Advice on using Arcade to make "mimic" the capabilities of Maplex in ArcGIS Desktop or is there a template or "standard" for label placement and scale ranges using Arcade?

To specify most of the labeling concerns road centerlines and their spacing, density and conflict with other feature layers in my Web map.


texts: Calculating distance to median point for a group of points in Model Builder In ArcGIS

texts:
I have two shapefiles. First of them, there is about 115 000 points (clients and their field size) All together there is about 16 226 different clients. In the other shapefile I have calculated the median point for all the clients. Now, I want to calculate the distance from each client field to its corresponding median point. How can i do it in Model Builder?? Doing it one by one takes a lot of time.

I really hope someone can help me.


texts:
I’m not an expert in GIS but recently I studied a bit, came on different sites. Among others QGIS, MangoMap, ArcGIS, … etc

A group of 6 small towns in Belgium wants to join efforts for bringing their graveyard online.

I’m wondering if the following procedure could be effective…

Take a photograph with coordinates from the sites (with a drone instead of an airplane…?)
Import the picture as a layer in AUTOCAD
In another layer draw Lines/Polylines over the existing  graves
Add drawing for known future implantation of graves (not much)
Export the AUTOCAD files to shape
Upload the files to f.e. MangoMaps  or another server
Place one POINT per deceased person on the graves (unique identification)
Make an Excel file with the same unique identifier per deceased person
Connect the Excel file with the files in MangoMap (upload)
Extra persons can easily be added in the file and re-uploaded

One town (mine) has a shapefile of the site. I can import this file in AUTOCAD and add the new graves. Some errors also have to be corrected.
Can also be done in QGis but not sure what’s easier.

Towns who have no shape files don’t need them. Just:

upload the picture
put the number of POINTS (deceased persons) on the graves (can be several points on one grave)
export the file as Excel (believe that’s possible)
add the deceased persons in the EXCEL file at the correct unique identifier
upload the Excel file

Could this be a good way to do it?


texts: Android Library to convert ESRI Shapefiles to GeoJSON

texts:
I'm building an Android App for display ESRI Shapefiles but unfortunately, there are only a few free SDKs which support Shapefiles.
Most of them working with GeoJSON, so now I'm looking for an Android Library to convert my Shapefiles to GeoJSON.

Is there any ready to use solution for this problem?
I don't need any standalone tool


texts: Figuring out what coordinates supposedly using 'NAD31 projection' refer to

texts:
I'm working with a GeoJSON file that contains sets of coordinates associated with street corners in the city of Montreal, Canada. The spec sheet of the API delivering those GeoJSON specifies that the coordinate system uses the NAD31 projection datum, which in spite of my search I doubt exist.

Because of the street names associated to the coordinates in the GeoJSON, I am able to deduce that the erroneous NAD31 coordinates 300247.509025, 5040897.253 refer to the WGS84 latitude/longitude 45.5079942,-73.5582652.

Does anybody know what might the geodetic datum NAD31 refer to, and how I can convert it to WGS84 coordinates?


texts: Interpolate points between coordinates for smooth animation in google maps or openmap

texts:
I have an application where we want to create an animated route based on recorded points.  I have an example that use the google maps driving directions, however, we already have points and they are not really on mappable roads.

I have animated the points I have already, but it's kind of choppy.  If I could interpolate some points between then I could create smoother animation.

I have location and time stamps for each so I could use that to determine how many points I want in between (for speed).

I am using google maps at the moment, but I don't really want to depend on too many specific google api calls because eventually I want to run this using something open source and stand alone (for demos off network).

Thanks


texts:
Looking for the equivalent of the Cut Fill Tool in QGIS. I have multiple DSM rasters from a drone that I would like to do a volume analysis on.

A cut-and-fill operation is a procedure in which the elevation of a landform surface is modified by the removal or addition of surface material. The Cut Fill tool summarizes the areas and volumes of change from a cut-and-fill operation. By taking surfaces of a given location at two different time periods, it identifies regions of surface material removal, surface material addition, and areas where the surface has not changed.

How can I do this in QGIS?


texts:
I have a selection of 30 features from a dataset in QGIS. I want to inspect every single feature one after another with the attributetable and the icon 'zoom to selected object'.
If I click on one row the selection is gone and only the clicked row is selected.
Is there any button which I have to press?
I don´t want to store the selection as a extra file.


texts:
In QGIS 3.4 I need a long string field to store one description (2000 characters) for each point. The project is based on *.shp files however they manage only 254.
I need to compile lots of formats in an atlas.
Any tips or suggestions for different file format?


texts: Apply formatting of symbols from ArcGIS Pro to ArcMap

texts:
I have a layer with about 11,000 points. Each point is tied to a point in another layer. The relationship is more or less a parent/children type, where all children tied to a parent is of the same color for their symbology. I am trying to change the outline for all 11,000 points without disrupting the color matching between children and parent points. When I go into the symbol properties and select "Properties for All Symbols", I can change the outline, but then it ends up changing the color of all points to the default monochrome.
An alternative that I've gotten to work is to do the same workflow in ArcGIS Pro by changing the border through the "Properties for All Symbols" which does not change the color and my points end up looking the way I intend. The issue is that I need to publish the service to ArcGIS for Server 10.5, which ArcGIS Pro can not do (only to Portal/AGOL). I was able to publish the two layers from Pro to my Portal and then add them to my ArcMap document as data from AGOL.
Visually, all looks good at the MXD level, but when I publish the MXD as a service to Server 10.5, the layers do not draw correctly, one layer will show, but the other will not. It seems like it's a layer order issue, because it's always the first ordered one that shows in the service and never the second one.
Does anyone know how I can solve this workflow issue or can show me a way that I can either change the symbology in ArcMap or export the layer properties in ArcGIS Pro to import into ArcMap?


texts:
I have written a code that is taking data from MODIS about NDVI, filter for March only , calculates the mean value for the whole month and clip it with geometry. The NDVI data is for period of years, for example, from 2000-2003.

The problem is that I need to run this process for few periods, 2000-2003, 2004-2007 and more and I dont want to use copy paste. So I want to create a loop but I'm coming from python and the exmples I have seen here in this forum confused me a little because I dont understand what the loop contains.

If someone can explain me the logic or how to put it inside a loop it will be very helpful.

This is the code I have worked on untill now:

var dataset = ee.ImageCollection('MODIS/006/MOD13A1')
.filter(ee.Filter.date('2000-03-01', '2003-03-31'))
.select('NDVI');

// filter the image collection to contain only images from March

var onlyMarch = dataset.filter(ee.Filter.calendarRange(3, 3,'month'));
print(onlyMarch);
// make a composite (mean image) of the images of March and clip to the geometry

var meanImage = onlyMarch.mean().clip(table).divide(10000);
print(meanImage);

var mean20002003 = meanImage.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: table,
  scale:500
});

//var ndvi = dataset.select('NDVI')
var ndviParams = {min: 0, max: 1, palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ]};

Map.centerObject(table, 7);
Map.addLayer(meanImage,ndviParams,'NDVI20002003');
print(mean20002003,'mean NDVI 2000-2003')


texts:
Can any one help me out to find route for three points. the points are start-via(Mid)-End. I want route to get create from start point to via(mid) point and via point to end point. Can any one can provide any help to me regarding this issue.
using this link:[http://workshop.pgrouting.org/index.html. How it will be done Sql query to retive same from php file in link above


texts: Iterate over every single value in a Raster GDAL C API

texts:
I am learning C, and I am trying to use the C API to iterate over the values of an array. So far I've been playing around printing the projection, Width and Height of the raster, but I would like to iterate over the values of a band I am reading as an array, here is my code (modified from enter link description herehttps://gdal.org/tutorials/raster_api_tut.html)

#include "gdal/gdal.h"
#include "gdal/cpl_conv.h"
int main()
{

    GDALDatasetH raster;
    GDALAllRegister();
    raster = GDALOpen( "/home/roger/Documents/98_GDAL_C++_snippets/data/Isle_wight.tif", GA_ReadOnly);
    if(raster==NULL)
    {
    printf("Invalid Raster Dataset");
    }
    else
    {
    GDALRasterBandH raster_band;
    int XSize, YSize;

    raster_band = GDALGetRasterBand( raster, 1);
    GDALGetBlockSize( raster_band, &XSize, &YSize);
    float *pafScanline;
    int nXSize = GDALGetRasterBandXSize (raster_band);
    int nYSize = GDALGetRasterBandYSize (raster_band);
    pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize);
    float band1[nXSize];
    band1 = GDALRasterIO(raster_band, 0, 0,GF_Read, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
}
}

In the last line: band1 = GDALRasterIO(raster_band, 0, 0, GF_Read, nXSize, 1, pafScanline, nXSIze, 1, GDT_Float32, 0, 0),

I receive the following message:
float band1[nXSize] : expression must be a modifiable value

I guess it is because I am trying to assing what comes from GDALRasterIO to an array, I thought that is what comes from GDALRasterIO, as the documentation states:

  Reading Raster Data

There are a few ways to read raster data, but the most common is via the GDALRasterBand::RasterIO() method. This method will automatically take care of data type conversion, up/down sampling and windowing. The following code will read the first scanline of data into a similarly sized buffer, converting it to floating point as part of the operation.

Is this first scanline of data an array I can iterate over? Or I am doing it in the wrong way?


texts: ArcMap Drawing Error One or More Layers Failed to Draw?

texts:
I am running ArcMap 10.3.1 on a 2014 MacBook Pro, Windows 10 via Bootcamp. I downloaded a cartographic shapefile off of the US Census Bureau: Urban Areas --- cb_2016_us_ua10_500k.zip. Adding the shapefile immediately resulted in the error:

  One or more layers failed to draw:

  cb_2016_us_ua10_500k:   There was a problem drawing the shapefile.
  Drawing aborted.

This is the first time adding shapefiles using ArcMap with this computer. Zooming in and out does not make a difference. I read this article and it does not seem to be the answer to my problem (https://support.esri.com/en/technical-article/000005082). I have not had issues adding geodatabases for a course I am taking. In the past I have had not an issue with shapefiles from the US Census Bureau. The shapefile does not appear to be too large either (attached is a picture of the folder).


texts: Problem with tool validation script in ArcGIS?

texts:
I have tried to write some validation code for a script in ArcGIS. This validation code more or less repeats what the script itself does, which is a simple cutting up of a long string by using the Python "split" string method, and then using each of the outcomes to define some parameters. This works great for most variables here, but for some reason the params 2 and 6 are being set to the full original string when I review the validation results in the ModelBuilder window by hovering over the variable. I have no idea why this is happening. Both of these parameters are defined as String output on the script. There is one more string output variable here (params 8), that displays a good value, and the only difference I see is the extra "replace". I have even tried to hard code the values for params 2 and 6, by setting it directly equivalent to a string, but to no avail, it still showed the full original unsplitted string.

What am I missing?

def updateParameters(self):
"""Modify the values and properties of parameters before internal
validation is performed.  This method is called whenever a parameter
has been changed."""
    self.params[1].value = self.params[0].value.split(';')[0]
    self.params[2].value = self.params[0].value.split(';')[1]
    self.params[3].value = self.params[0].value.split(';')[2]
    self.params[4].value = self.params[0].value.split(';')[3]
    isExportFeatureClassStr = self.params[0].value.split(';')[4]

    if isExportFeatureClassStr == "True":
        isExportFeatureClass = True
        flowControlBoolean = False
    else:
        isExportFeatureClass = False
        flowControlBoolean = True

    self.params[5].value = isExportFeatureClass
    self.params[6].value = self.params[0].value.split(';')[5]
    self.params[7].value = flowControlBoolean
    self.params[8].value = self.params[0].value.split(';')[6].replace(",",";")

return


texts:
With Python 3 and GDAL package, I am trying to merge two shapefiles that come with databases.

They come with the following files:

grade_id14.cpg
grade_id14.dbf
grade_id14.prj
grade_id14.sbn
grade_id14.sbx
grade_id14.shp

My code goes as it follows :

from osgeo import gdal, ogr, osr

test = ogr.Open('grade_id13.shp')
test2 = ogr.Open('grade_id14.shp')

lyr1 = test.GetLayer()
lyr2 = test2.GetLayer()
unionc = ogr.Geometry(ogr.wkbMultiPolygon)
print(lyr1.GetGeomType()) #it shows 3

for feat in range(lyr1.GetFeatureCount()):
    fit1= lyr1.GetFeature(feat)
    geom1= fit1.GetGeometryRef()
    unionc.AddGeometry(geom1)

for feat2 in range(lyr2.GetFeatureCount()):
    fit2= lyr2.GetFeature(feat)
    geom2= fit2.GetGeometryRef()
    unionc.AddGeometry(geom2)

union= unionc.UnionCascaded()

unionTest = geom2.Union(geom1)

# Create datasource.
driver = ogr.GetDriverByName('ESRI Shapefile')
try:
    ds_out = driver.CreateDataSource('test.shp')
except:
    raise Exception("File already exists.")

# I believe here is the problem, when I try to get the layers from the united polygon,
# doesnt work
lyrU = union.GetLayer()

srs = lyrU.GetSpatialRef()
geom_type = lyrU.GetGeomType()

# Create empty layer.
layer_out = ds_out.CreateLayer(name = 'testLayer', srs = srs, geom_type = geom_type)

# Copy fields <-- not required.
layer_defn = lyrU.GetLayerDefn()
n_fields = layer_defn.GetFieldCount()
for i in range(n_fields):
    layer_out.CreateField(layer_defn.GetFieldDefn(i))

My ultimate objective is to unify several shapefiles.


texts: Append object from another layer to the new layer through form (drop-down list)

texts:
I'm trying to get to know is it possible to append object from another layer to the new layer through form (drop-down list)? Do I need to create an application? If yes please tell me what language would be the best for this kind of problem?

Example: I created new layer called Companies_objects. I have another layer (all_objects) with a lot of objects in town and some of theme belong to specific companies. I would like to append some Companies that interesting me (by it's names) from the drop-down list from all_objects layer when I'm editing new Companies layer. Is it possible?


texts: Force geometry type in mixed geometry WFS layer in QGIS?

texts:
I have an external WFS source serving mixed geometry types - in other words, I have no control of the source. When adding the layer QGIS sets the geometry type to line and thus I am unable to style the points and polygons also contained within the layer.

Nathan W states in a comment that a layer in QGIS can contain only one geometry type. If I was able to explicitly set the geometry type on a layer, I could add the layer three times and use one layer for each geometry type.

Is it possible to force the geometry type on a (WFS) layer in QGIS?

Or is there another method I could use to correctly style all the features in my WFS layer?


texts: IFeature.set_Shape The coordinates or measures are out of bounds

texts:
I'm trying to add a point to my map, but something with the way I used shp2sde create caused me to have a incorrect spatial extent or something. So I'm getting "out of bounds" error.

How do I go about fixing this in ArcView 9.3.1? (I'll have to use command prompts probably or recreate).

My Map Server Info:

Intial Extent:

XMin: -197.114430117841
YMin: -95.6461199683022
XMax: 197.694572740662
YMax: 184.916970645712
Spatial Reference: 4326

Full Extent:

XMin: -198
YMin: -99
XMax: 198
YMax: 99
Spatial Reference: 4326

My Layer info:

Extent:

XMin: -180
YMin: -90
XMax: 180
YMax: 90
Spatial Reference: 4326

Here's my log file from what I'm trying to INSERT into the feature layer:

Check geometry point type=esriGeometryPoint
about to create feature
Create shape(X=38.9300285250005,Y=-74.9343213739996)
Com_point_spatialreference=GCS_WGS_1984
at ESRI.ArcGIS.Geodatabase.IFeature.set_Shape(IGeometry Shape)
 The coordinates or measures are out of bounds.

My sdelayer -o describe_long:

Spatial Column .......: SHAPE
Layer Id .............: 21
SRID .................: 5
Minimum Shape Id .....: 1
Offset ...............:
  falsex:          0.000000
  falsey:          0.000000
System Units .........:       1000.000000
Z Offset..............:          0.000000
Z Units ..............:          1.000000
Measure Offset .......: <None>
Measure Units ........: <None>
XY Cluster Tolerance .:          0.002
Spatial Index ........:
  parameter:    SPIDX_GRID,GRID0=1,FULL
  exist:        Yes
  array form:   1,0,0
Layer Envelope .......:
  minx:      -180.00000,        miny:       -90.00000
  maxx:       180.00000,        maxy:        90.00000
Entities .............: np
Layer Type ...........: SDE
Creation Date ........: 03/21/12 11:04:05
I/O Mode .............: LOAD
Autolocking ..........: Enabled
Precision.............: Basic
User Privileges ......: SELECT, UPDATE, INSERT, DELETE
Coordinate System ....: GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1
984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292
5199433]]

Layer Configuration ..: DEFAULTS


texts: Multivalue processing in the context of the modelbuilder

texts:
I am new to python so please bear with me. I have the following script that is supposed to get a multivalue parameter of several shapefiles as an input and then run the unionanalysis on all these shapefules to create a single feature class or shapefile output. The script seems to run fine however only processes up to step 2 of 16 and does not create any output file. I do have the output file specified however and all the parameters are in the right order with the input accepting multivalue variables and both having feature class as their data type. I need to precise that this is within the ModelBuilder framework...

Here is the script:

    import arcpy

    myList = arcpy.GetParameterAsText(0)

    unionResult = arcpy.GetParameterAsText(1)
    arcpy.Union_analysis(myList, unionResult, "ALL", "", True)

To further clarify:

Parameter (0) is an input, feature class, multivalue (yes)
Parameter (1) is an output, feature class, multivalue (No)


texts: Can I tell QGIS to load/show only a part of an image?

texts: can I tell QGIS to load/show only a part of an image without to extract this region as an extra/ new image?


texts: I am customizing the opentripplanner webapp according to my needs. it has 5 different base layers. but in code (in map.js and OpenLayersUtils) i can't find where it adds these layers to map.


texts: Is there a way to disable the FieldCalcBar thingy, which appears in attribute table when the layer is in editing mode? I suppose there should be a checkbox in Settings -> Customization... dialog window, but so far haven't found the right one.


texts: Is there a way to export a point layer from QGIS to a csv file?

texts: I have saved waypoints from my Garmin as gpx and have successfully imported them into QGIS. Is there a way to now export the layers into a csv file with coordinates? I've tried the export to csv but that does not give me the coordinates.


texts: Editor toolbar disappears whenever mxd is re-opened

texts:
We just switched to ArcGIS 10.3.1 from 10.2.2 at work. I've noticed that when I start to edit a feature, the Editor toolbar appears. I dock, customize it a bit, make and save some edits, and close the mxd. When I open the map again, the Editor toolbar is gone (Advanced Editing toolbar that I added remains). I start editing again, Editor toolbar pops up, minus the customization.

I can't find a setting to keep it as I added it. Never seen this behavior in previous versions of Arc. Is there a way to make it behave like previous versions?


texts: Temporarily hiding features that are not selected at the moment in Leaflet

texts:
How can I set focus on the selected feature by additionally and temporarily hiding other features that are not selected at the moment.
Here is my code:
var highlight = {
       'color': '#149c14',
       'weight': 4,
       'opacity': 1
   };

$.getJSON("/data/geojson/gen.json", function(data){
           var genlayer = L.geoJson(data,{
           style:function(feature){
               return {color:"red"};
           },
           onEachFeature:function(feature, layer){
               layer.on('click', function(e){
                    //do functions to selected feature
                    selectedFeature = e.target;
                    layer.setStyle(highlight);
               });
           }
      });
    roads.addLayer(genLayer);
   });


texts: Is the TIN created out of several cross sections different from merging the TINs created out of each set of two cross sections?

texts:
I have a river with cross sections that have the water surface elevation values. I am trying to create TIN from the water surface values and then compare it with the topo (to see the flooded areas). I wonder what would be the difference between the TIN created out of several cross sections is different from merging the TINs created out of each set of two cross sections? Does the TIN creation process includes any sort of generalization/ approximation/ getting average when it comes to several features or it is dependent on the adjacent features only?
My tests on the river shows no distinct difference between the resulting floodplain but I am not quite sure it is the same everywhere. Any thought?


texts: Shapefile with overlapping polygons: calculate average values

texts: I have a very big polygon shapefile with hundreds of features, often overlapping each other. Each of these features has a value stored in the attribute table. I simply need to calculate the average values in the areas where they overlap. I can imagine that this task requires several intricate steps: I was wondering if there is a straightforward methodology. I’m open to every kind of suggestion, I can use ArcMap, QGis, arcpy scripts, PostGis, GDAL… I just need ideas. Thanks!


texts: Identifying LiDAR tiles that are located within boundaries of AOI

texts:
I am attempting to replicate the process of creating digital elevation models from raw LiDAR data and have been provided with LAS files covering an approximate 530 square mile county area. However, I only need to work with about 1/4th of this coverage, which is still at least 30+ las files.

What is the most efficient way to identify the files that are located within the boundary of my AOI (area of interest)?

The company that delivered the LiDAR data did include an ArcMap map document that contains a point cloud LAS dataset file but I'm not sure how/if I can manipulate that along with my boundary polygon.


texts: Extending (buffering) set of polygons into empty space only using Shapely?

texts:
I have a set of polygons of which I would like to extend the boundaries a bit into empty space, using Shapely. They are not allowed to overlap each other at all. For example creating buffers and then removing the overlapping areas is not an option, since this will remove a part of the original areas.

The following image (partially) indicates how I would like the result to be like:

I have retrieved the data through a WFS service, and parsed it into Shapely polygons. I've also created buffers. Like so, in a loop:

# in this case, the feature's geometry only has an exterior boundary
if len(bag_polygon) == 1:
    bag_polygons.append(Polygon(bag_polygon[0]))
    bag_polygons_buffered.append(Polygon(bag_polygon[0]).buffer(1, cap_style=3, join_style=2))

# in this case, it has interior boundaries as well
elif len(bag_polygon) > 1:
    bag_polygons.append(Polygon(bag_polygon[0], bag_polygon[1:]))
    bag_polygons_buffered.append(Polygon(bag_polygon[0], bag_polygon[1:]).buffer(1, cap_style=3, join_style=2))

Putting them in a .shp with Fiona, I have added a vital attribute list to the geometries.

shp_schema = {
    'geometry': 'Polygon',
    'properties': {'bag_id': 'int'},
}

with fiona.open("bag_subset.shp", "w", crs=from_epsg(28992), driver="ESRI Shapefile", schema=shp_schema) as bag_shp:
    for i in range(len(bag_polygons)):
        bag_shp.write({
                "geometry":mapping(bag_polygons[i]),
                "properties": {"bag_id":bag_ids[i]}
                })

with fiona.open("bag_buffered_subset.shp", "w", crs=from_epsg(28992), driver="ESRI Shapefile", schema=shp_schema) as bag_shp:
    for i in range(len(bag_polygons)):
        bag_shp.write({
                "geometry":mapping(bag_polygons_buffered[i]),
                "properties": {"bag_id":bag_ids[i]}
                })

Lastly, I've thought about taking the difference between the shapefile that includes the buffers and the shapefile with the original geometries, to subsequently join this difference with the original geometries. However, there is still the problem that the buffer zones overlap themselves.


texts: Creating setback buffer for lines using QGIS?

texts:
I've attached a little diagram showing what I want to do. Using the directional buffer tool results in the bottom left output, which is not what I want, which is the center one.

Regular buffer tools result in the diagram on the right, which would have been one step towards the solution I'd have gone for: using the coastlines shapefile as the other layer in an intersection operation that would result in what I want.


texts: How to manipulate SDF data via java/jsp (MapGuide open source Web API)?

texts: I am relatively new to MapGuide open source technology. I would like to parse and modify data that is in SDF format so that my map displays this updated feature data (dynamically). How can i do so? I have MapGuide installed and configured on Linux. I went through the MapGuide Web API documentation but there were no examples on how to manipulate SDF data. Currently i am generating the SDF data in a different environment (in a Windows OS domain) by using scripts to query a spatial database and the query generates SHP files. I then transform those SHP files into SDF files (using scripts) before transferring them to the MapGuide server repository in the Linux OS environment. On the Windows OS end, I am able to edit the SHP files (in particular the dbf file just using a spread sheet) to change the values for the particular feature source. However this process is static and so i was hoping to be able to do so programmatically (using java / jsp) by parsing the SDF files for the particular feature source in the Linux domain. Thanks


texts: GeoTIFF appears only as a black solid image

texts: I need to open a GeoTIFF file (high resolution aerial image) but for some reason it appears only as a black solid image and nothing else. I've tried opening with QGIS and Windows photo viewer so far.  What can I do to solve this?


texts: Cannot get number of max pixels in Google Earth Engine

texts:
I have a code that calculates the max number of pixels among all the rasters it retrieves and then creates new imagecollection with only rasters that contains  rasters in the size I want. For some reason, , today I have tried to do this on different dataset- sentinel 2, but it didn't work  and I get null as number of pixels.

I suspect that it's because on the prexious dataset I have used dataset that had only one band image- NDVI, and now I have more, but I don't really know and without this number I can't run the whole script.

This is the code I have:

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

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

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

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

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

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

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

// Get the number of images.
var count = clippedCol.size();
print('Count: ',count);

// //Filter according to number of pixels

var ClippedCount = clippedCol.map(function(image){
  var countpixels = ee.Number(image.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: geometry,
  crs: 'EPSG:4326',
  scale: 30,
  //maxPixels:21176024,
  }));

  return image.set('count',ClippedCount);
});

print('Count:',ClippedCount);

//Here I get null in print
var max = ClippedCount.reduceColumns(ee.Reducer.max(),  ["count"]);
print(max.get('max'));
var max_number=max.getNumber('max');

I need to filter the images in the end according to this number, I am not sure how to solve it .

*I can't upload the original geometry but it didn't work on any geometry I have tried so you can try it on whatever polygon you draw


texts: GE Smallworld PNI 4.3 forgot Administrator Password

texts: I forgot the Administrator password for my company's Smallworld PNI installation, we had been making due with just the user accounts but now the user accounts don't log in anymore, is there anyway to reset the Admin password


texts: Building OpenLayers map with many layers efficiently?

texts:
I have built a WMS parser which I am using to retrieve all layers from a WMS url to be able to automatically build a OL map with those layers inside.

Actually, I see I can do it in two modes:

MODE 1: making a ol.source.ImageWMS instance for each layer and append it to a layer_list array (to be used as ol.Map's params parameter)
MODE 2: making only one instance of ol.source.ImageWMS with all the layer names inside the params:LAYERS parameter.

To my understanding, the MODE 1 should make a new request for each layer to the same WMS url, and this have the disadvantage of slowing the process, but the advantage of maybe handling the layers easier (e.g turing them on and off, zooming to one specific layer, etc.)

MODE 2 however migt be harder to handle in terms of single layer, but it's only one WMS request.

I am new to OL.

Which mode is more efficient?


texts: What is a convenient way to write OSM data to a .shp file in python

texts:
I wish grab OSM data that is within a specified bounding box and matches specific key/value combinations. I then want to write this to a .shp file for further processing using, probably, the QGIS api.

I have successfully acquired OSM data within the bounding box using overpy and have it within an overpy.Result object. I just can't find a way to write it to a .shp file?

Can anyone direct me to a suitable fragment, preferably one that uses fiona?


texts: How can I rotate a raster layer (global land-mass) to make the 180 deg longitude the center of the map (using QGIS)? 'rotate' function in R is simple but not accurate. Raster has 'longlat' projection and WGS84.


texts: QGIS 3.10 export layer to KML reverses coordinates

texts:
I'm using QGIS 3.10. I generated a contour map from a DEM. It looked good. Now, I wanted to check the contours in Google Earth. From the contour layer in QGIS, I selected Export - Save Features As - KML. I ran the export and a KML layer was generated as a file. Problem is, Google Earth expects the coordinates to be longitude, latitude but the exported file has them as latitude, longitude. When the KML file is displayed in GE, it zooms to Antarctica.

The schema used by QGIS is KML 2.2 which suggests the required coordinate sequence is longitude, latitude. But QGIS is sequencing latitude, longitude.

Is there a way in QGIS to specify the sequence to be longitude, latitude?

Excerpt from the KML file generated by QGIS:

<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document id="root_doc">
<Schema name="Watersprite_contours" id="Watersprite_contours">
    <SimpleField name="fid" type="float"></SimpleField>
    <SimpleField name="ID" type="int"></SimpleField>
    <SimpleField name="ELEV" type="float"></SimpleField>
</Schema>
<Folder><name>Watersprite_contours</name>
  <Placemark>
    <Style><LineStyle><color>ff0000ff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>
    <ExtendedData><SchemaData schemaUrl="#Watersprite_contours">
        <SimpleData name="fid">1</SimpleData>
        <SimpleData name="ID">0</SimpleData>
        <SimpleData name="ELEV">800</SimpleData>
    </SchemaData></ExtendedData>
        <LineString><coordinates>50.0173281602732,-122.806140688556 50.0172702778432,-122.806140688556 50.0172606307647,-122.806050857117 50.0172702778432,-122.80587119424 50.0173281602732,-122.80587119424</coordinates></LineString>
  </Placemark>


texts: How to randomize colors for large datasets in QGIS new symbology?

texts:
Someone has already asked how to randomize colors by class in QGIS new symbology, but the answer there does not work with datasets containing more than 100 unique values since the Random option only allows you to create only 100 colors. Using that technique, the features get categorized using the equal breaks method making the color the same for many consecutive features instead of assigning one per feature and then starting over with the 100 colors again after they have all been used once. My goal is to differentiate features in a polygon layer with many similarly located and named features, e.g., Woodlands 1 and Woodlands 2 are next to each other and so need to be different colors.

Is there a way to make sure that each consecutive feature is a different color from the immediately preceding and following feature? In the old symbology, you can just use "Randomize" to achieve this. Is there a way to do it in the new symbology? Thanks for any ideas or help!


texts: How are EGM96 and WGS84 related to each other?

texts:
As far as I understand the EGM96 defines the Geoid, where as the WGS84 Standard defines the Ellipsoid.

Is the ellipsoid defined in the WGS84 standard defined in a way to maximize the congruency with the the geoid defined in the EGM96 standard?


texts: Overlap area of features from same feature class

texts: Is there a way to calculate the area of overlap between features (polygons) from the same feature class. See figure below as illustration for example. There are 3 features with an ID and area and I want to calculate the area where 1 overlaps with 2, 1 overlaps with 3, and 2 overlaps with 3. The idea is to get a sort of cross-validation table with polygons 1, 2 and 3 across each row and column with the diagonal as 1 (i.e. 100%) and the off diagonal as the proportion overlap with other polygons. Keeping in mind that this is all one feature class or shapefile. I'm thinking the Identity tool in ArcGIS ...


texts: How can I turn islands into individual features in QGIS?

texts:
I have what seems like a pretty straightforward task in QGIS: I have a shapefile that, in some cases, groups multiple separate landmasses into a single feature.  For example, the Hawaiian islands are all one feature.  How can I split those up into multiple features?  This seems like it should be easy but it's proving impossible.

Using the "Split Feature" tool doesn't do anything when I try to draw a split line between the islands, or surround one of them.  The "Singlepart to Multipart" option in the vector tools might do the job, but it's not selective, it just splits up EVERY multipart feature, when I only want to split up a few of them.

How can I pick a polygon that's one part of a multipart feature, like a single island in a chain, and split that out as a separate feature?


texts: Python: How to specify crs of raster file by Python/rasterio?

texts:
I use the rasterio to cut out roofs of buildings, roofs are polygons. I would like to set up crs of felled roofs. Is it possible to do that by rasterio? I can set up crs of layer, layer.crs = 'epsg:2180', but how can I do that for raster file?
My code to cutting out roofs of buildings:
def cutting_out(df, start, end):
    for index, row in df.iloc[start:end].iterrows():
        if len(row['1_5000_pi025']) == 1:
        #1url
            coo = shapely.wkt.loads(row['geometry'])
            multiplier = int(row['liczba_kon'])
            coo_buffer = coo.buffer(3 * multiplier)
            lokalnyid = row['lokalnyid']
            tiff = str(row['1_5000_pi025'][0])
            with rasterio.open(f'/mnt/qnap/geo_tiff_1_5000_p025/{tiff}') as src:
                out_image, out_transform = rasterio.mask.mask(src, [coo_buffer], crop=True)
                out_meta = src.meta
                out_meta.update(
                    {
                        "driver": "GTiff",
                        "height": out_image.shape[1],
                        "width": out_image.shape[2],
                        "transform": out_transform,
                    }
                )
                with rasterio.open(f'/mnt/qnap/geo_tiff_roofs_testy3/{lokalnyid}'+'_025'+'.tif', "w", **out_meta) as dest:
                    dest.write(out_image)


texts: ArcGIS Server - Get map of all folders, services, layers etc

texts:
Is there a way to query a REST service on an ArcGIS server and get JSON object of all folders, services, layers and other resources.

I am writing my own queries but it is becoming rather cumbersome and I am wondering if something is already built in

var servicesRoot='https://my.gisserver.com/gisServer/rest/services';

$.ajax({
  url: servicesRoot+'?f=pjson',
  dataType: "jsonp",
  success: function(response) {
    //this logs all the folder in this service
    console.log(response.folders)
  }
});


texts: Empty result when writing two band TIFF using GDAL/Python

texts:
I am trying to write two band TIFF using GDAL/Python. I get empty result although the I see the array should be written in my variables.
    driver = gdal.GetDriverByName("GTiff") # make driver with the geotif format

outdata = driver.Create(outFileName, rows, cols, bands, gdal.GDT_Float32) # parameters: output name, rows, col, bands, data type
outdata.SetGeoTransform(a.GetGeoTransform())##sets same geotransform as input
outdata.SetProjection(a.GetProjection())##sets same projection as input
for i in [0,1]:
    outdata.GetRasterBand(i+1).WriteArray(arr_out[:,:,i])
    outdata.FlushCache() ##saves to disk!!


texts: Rendering gridded time-series data with Leaflet.js

texts:
Background

My colleague and I have developed a web-app for visualizing gridded time-series data (in the arctic) using leaflet.js with webgl-heatmap.

The source data is in netCDF format, which we have been processing to geojson for each parameter attribute and and time-step. This data is then accessed via an api depending on what the user has selected.

Problem

It is difficult to stylize a color gradient for the heatmap to accurately reflect levels of the data, and WebGL presents problems for various browsers (and of various versions).

Inspired by the geojson-vt plugin, our idea now is to render each cell as a geojson polygon, and update the style (i.e. color and opacity) of those cells according to their values for the current data selection—possibly processing this as vector tiles to improve perfomance.

Proposed work-flow

netcdf → PostGIS
grid cell geometries → vector tiles
user data selection & time-step → query PostGIS for values → style grid cells

Question

What workflow and tools would be best suited to achieve this?


texts: Locating where two lines of defined length meet, maybe with QGIS Plug-in

texts:
On a site, we recorded various pits by measuring the two distances to two known physical objects (e.g. 20 metres to pit from physical object A, and 10 metres to pit from physical object B).
I am now adding the location of these pits to a polygon layer I'm working on but it is cumbersome measuring 20 metres from object A roughly, then 10 metres from object B roughly, repeating it a few times to get exactly where the two lengths meet up.
Is there a method or a plug-in that would do this easily? (There can be only two possible points where they can meet, of course.)


texts:
My PostGIS database contains OSM data from my country. For transportation reasons I want to export a spatial view I created from that data to a SpatiaLite database.

Is there any way to directly extract my PostGIS view to a SpatiaLite db without resorting to a man-in-the-middle (shapefile then SQLite) solution?

Im using Ubuntu Linux, with PostGIS 2.1 GDAL 191.

EDIT:

Here's a single cmd to export to shapefile:

ogr2ogr -f "ESRI Shapefile" <name-of-folder-to-be-created> \
PG:"dbname=dbname user=username password=pass" \
-sql "<sql query eg> select * from database where <whatever>"

If you get a getsrid(geom) error, load in your database the "legacy.sql"

psql -f /path/to/legacy.sql


texts: Limiting feature class parameter to specific geodatabase?

texts:
I have a tool that will update the geometry of a feature class based on the geometry of another, using a common identifier.

It works great but the way I have it set up, you have to first select the gdb holding the two feature classes, then select each feature class.

Is there a way, once the gdb (which defines the workspace) is selected, to then limit the options for the two feature class parameters to those found in the selected gdb?

When I run it, I know they have to be in the same gdb, but I'd like to share it with colleagues as a tool and I'd like it to be at least a little fool-proof, considering it was fool-made.

Here's the code:

import arcpy

arcpy.env.workspace = arcpy.GetParameterAsText(0)
UpdateData = arcpy.GetParameterAsText(1)
UpdateField = arcpy.GetParameterAsText(2)
ProjectData = arcpy.GetParameterAsText(3)
ProjectField = arcpy.GetParameterAsText(4)

geometries = {key:value for (key,value) in arcpy.da.SearchCursor(UpdateData, [UpdateField, 'SHAPE@'])}

notfound = []

with arcpy.da.UpdateCursor(ProjectData, [ProjectField, 'SHAPE@']) as cursor:
    for row in cursor:
        try:
            row[1] = geometries[row[0]]
            cursor.updateRow(row)
        except:
            notfound.append(row[0])


texts: How to reproject all geometries columns in a PostGIS schema

texts:
Is there's a way to reproject all geometries columns of a postgis schema given the same source SRID to the same destination SRID?
Exemple :
schema.table1 has 2 geometries columns geom and geom2 with the same SRID 900913
schema.table2 has 2 geometries columns geom and geom2 with the same SRID 900913
etc..
I want to project all those columns in SRID 2154.
I tought there's could be a simple way by using the public.geometry_columns view but I don't know how.


texts:
I want to create a buffer on a certain layer that will have radios of n meters.

I go to vector->geoprocessing tools -> and then nothing buffer(s) is missing I have enabled the Spatial Query Plugin

Any idea how to fix this?


texts: Label only the first letter of string value

texts:
I have a field called 'Genus' where there are string values, for example word 'Pinus'.
I want my label to show only the first letter of this word, so in my example 'P'.
I was trying to use Function Editor. I was contriving with "field.row[0]", but I'm not sure with a rest of syntax.
Maybe there is another solution for it?


texts: Creating hiking track based on Open Cycle Map (OCM) data?

texts:
I need help.
In QGIS i’m trying to create GPX hiking route having Open Cycle Map as background layer. Like in brouter selecting points on map and to get route following routes on OSM.
Successfully imported OSM into QGIS but sometimes when i select parts of track it selects all branch and I need just part of that branch.