texts: Summarizing flatness of county with statistic Elevation variance per square mile

texts:
I need to classify every county in the US as either hilly, moderately hilly, or flat. Accuracy is not critical, just a rough classification. I used R and followed the advice in this answer to import elevation and counties and calculated the zonal variance. This was inherently biased towards larger counties since they have more surface for hills, so I need to divide by area.
If I divide variance by square miles does that reasonably represent hilliness?
The preliminary results seem decent, with the highest values being places like Colorado, Hawaii, Montana, etc.


texts: PostGIS do not want GeometryCollection in output

texts:
I am intersecting two layers. 1. truckmun which has multipolygon and polygon geometries in it 2. streamsw which is just multistring. The SRID is 3424 for both layers.

Here is my intersecting query

create table streamsw2 as
select truckmun.mun as mun,truckmun.county as county, st_intersection(truckmun.geom,streamsw.geom2d) as geom
from truckmun,streamsw
where st_intersects(truckmun.geom,streamsw.geom2d);

the query runs but spits out geometry types of St_geometrycollection, st_linestring and st_multilinestring

to get an idea of the layers being intersected the truckmun layer was created by buffering a network and using the st_intersection function with municipaliites. streamsw is just a multilinestring streams layer that at certain locations spatially intersects the truckmun layer.

My goal here is to 1. understand why the geometry collection exists 2. not have the geometry collection in the new table of streamsw2


texts:
I want to add Yandex satellite imagery map as layer by using WMS server. But I cannot find one. The only link I could find is TMS server (which is working fine): https://sat01.maps.yandex.net/tiles?l=sat&v=3.379.0&x={x}&y={y}&z={z}

My questions are:

Is there any WMS servers with Yandex maps?
Is there any workaround to use TMS server as WMS?


texts: Implementing 'Add feature' functionality programmatically to QGIS?

texts:
I have a script that is adapted from the book
Menke - Mastering QGIS and used to draw a polygon directly on the map canvas.
Now I want to expand upon that script and mirror the 'add feature' functionality. I.e. the script is supposed to take the coordinates that are stored in a list, transform them into a geometry and add the feature.

I tried the following and while I'm not getting any errors, I continue to have only a single rubber band connecting all my polygons.

Here's the code:

from PyQt4.QtCore import *
from PyQt4.QtGui import QColor
from qgis.core import QGis, QgsMapLayerRegistry
from qgis.gui import QgsMapToolEmitPoint, QgsRubberBand, QgsMapTool

# Create new virtual layer
vlyr = QgsVectorLayer("Polygon", "temporary_polygons", "memory")
dprov = vlyr.dataProvider()

# Add field to virtual layer
dprov.addAttributes([QgsField("name", QVariant.String),
                     QgsField("size", QVariant.Double)])

vlyr.updateFields()

# Access MapTool
previousMapTool = iface.mapCanvas().mapTool()
myMapTool = QgsMapToolEmitPoint( iface.mapCanvas() )

# create the polygon rubber band associated to the current canvas
myRubberBand = QgsRubberBand( iface.mapCanvas(), QGis.Polygon )

# set rubber band style
color = QColor(78, 97, 114)
color.setAlpha(190)
myRubberBand.setColor(color)

# create empty list to store coordinates
coordinates = []

# Create Feature
poly = QgsFeature()

# define function draw rubberband
def drawBand( currentPos, clickedButton ):
    iface.mapCanvas().xyCoordinates.connect( drawBand )

    if myRubberBand and myRubberBand.numberOfVertices():
        myRubberBand.removeLastPoint()
        myRubberBand.addPoint( currentPos )

    if clickedButton == Qt.RightButton:
        pass

# define function mouseclick
def mouseClick( currentPos, clickedButton ):
    if clickedButton == Qt.LeftButton:
        myRubberBand.addPoint( currentPos )
        coordinates.append(currentPos)
        print coordinates
    if clickedButton == Qt.RightButton:

        dprov.addFeatures([poly])
        poly.setGeometry(QgsGeometry.fromPolygon ([coordinates]))
        vlyr.updateExtents()
        iface.mapCanvas().refresh()
        QgsMapLayerRegistry.instance().addMapLayers([vlyr])
        del coordinates[:]

myMapTool.canvasClicked.connect( mouseClick )
iface.mapCanvas().setMapTool( myMapTool )

Also, I tried to have a look at the code of the qgis 'add feature' tool, but was unable to find it within the repository. Perhaps someone could point me into the right direction here as well?


texts: Not able to close custom edit form in ArcPad 10.2.2 using ok button

texts:
Summary of problem: I have problem with ArcPad in that the "ok" button on my custom edit forms is unresponsive when I am trying to add a new feature.

Project set-up: I am using ArcMap 9.3.1 and ArcPad 10.2.2. I have three point feature layers in a geodatabase within ArcMap ("centres","probes" and "active") which I am checking out in .axf format for use with ArcPad.  Each of the three layers has an edit form in .apl format which I am linking to when I check the data out. I have a small piece of vbscript code in a .vbs file for doing some additions and subtractions to populate some fields based on other fields.This is called by both "centres" and "probes". The .apl files and the .vbs file appear to be correctly incorporated into the .axf file when I am checking data out. I also have a dropdown list in .dbf format which is accessed by the "centres" layer; I have been adding this manually to the relevant folder on my handheld as I could not see a way for this to be incorporated within the .axf. I also have a fourth point feature layer ("waypoints") which I am adding to my ArcPad maps as a separate .shp shapefile (because .axf format cannot cope with the autoincrement option which I require).

Layer editing: Both the "centres" and "probes" layers have predetermined points whose attributes are amended or entered in the field via the edit form. I also need the option to add new point features to each of these layers in the field and to enter attributes for these features. The "active" layer and the "waypoints" layers have no predetermined points, all features are added in the field.

Previous problems: I previously had several serious problems with the .axf file when using ArcPad 8 and ArcPad 10. These included all symbology for a layer disappearing when editing a feature, error messages ("Unable to save data due to error 80040E014") and the spontaneous, error-message-free loss of attribute data from a seemingly arbitrary number of previously visited point features. This latter (major) issue appeared to be caused when rapidly changing the selected feature. These problems all seem to have been fixed by updating to ArcPad 10.2.2, however this update has created a new problem which did not occur when using the previous versions.

Current problem: When creating a new feature in the "centres" or "active" layers on a handheld, the edit form correctly opens and can be filled in. However, clicking the "ok" button fails to close the form and save the new feature. The button flashes when clicked but the forms fails to close. I can close the form by clicking the red X/cancel button but not even the geometry of a new feature is saved. Furthermore, if I then open the edit form for an existing point in that layer, I find that the "ok" button has also ceased to function there. Closing and re-opening the map restores the function of the "ok" button for existing points. Trying to create new features within these layers in ArcPad on the PC simply crashes the application.

NB: I do not have any problems adding new features normally to the "waypoints" .shp file and I can add points to the "probes" layer which is in the .axf file, although this success seems inconsistent.

Solutions? The ability to add features in the field is fairly fundamental to ArcPad as is the ability to use .axf files, so I would rather not have to move the whole project over to standard shapefiles as this would prevent the easy collation to data from multiple fieldworkers.


texts: Making base map from image with known corner coordinates?

texts:
I'm trying to make a map in Openlayers using an image of a crater from mars (specifically this image) so I can add markers/lines/etc on it.

I know the min/max latitude/longitude for this image, so how would i go about turning this into a map i can zoom/pan/add markers/etc to?

I've used MapTiler but the free version has watermarks, the coordinates seem to be off and it's got the earth as it's basemap:


texts: Heatmap not working on FeatureLayer (ArcGIS JS API 4.8)

texts:
I tried my best but didn't get success to show heatmap on FeatureLayer with latest ArcGIS JavaScript API 4.8. My Code is located here.
I used HeatMapRenderer to render the point on the map as the heat map. But somehow it is not working at all. I think I am missing something which would be an easy fix.


texts: Create a bounding box from point shapefile

texts:
I have a BigTiff image and an associate point shapefile with features. I would like to know how can I create a bounding box from that shapefile to create/clip images around that points to train a model.

EDIT:
I'm trying to do it automatically in a Jupyter Notebook and I've followed this threat but all exported images are completely black.

from osgeo import ogr, gdal

InputImage = 'XXX.tif'
Shapefile = 'XXX.shp'
RasterFormat = 'GTiff'
PixelRes = 0.5
VectorFormat = 'ESRI Shapefile'

# Open datasets
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()

VectorDriver = ogr.GetDriverByName(VectorFormat)
VectorDataset = VectorDriver.Open(Shapefile, 0) # 0=Read-only, 1=Read-Write
layer = VectorDataset.GetLayer()
FeatureCount = layer.GetFeatureCount()
print("Feature Count:",FeatureCount)

# Iterate through the shapefile features
Count = 0
for feature in layer:
    Count += 1
    print("Processing feature "+str(Count)+" of "+str(FeatureCount)+"...")

    geom = feature.GetGeometryRef()
    minX, maxX, minY, maxY = geom.GetEnvelope()
    xmin, ymin = minX - 50, minY - 50
    xmax, ymax = maxX + 50, maxY + 50

    # Create raster
    OutTileName = str(Count)+'.SomeTileName.tif'
    OutTile = gdal.Warp(OutTileName, Raster, format=RasterFormat, outputBounds=[xmin, xmax, ymin, ymax], xRes=PixelRes, yRes=PixelRes, dstSRS=Projection, resampleAlg=gdal.GRA_NearestNeighbour, options=['COMPRESS=DEFLATE'])
    OutTile = None # Close dataset

# Close datasets
Raster = None
VectorDataset.Destroy()
print("Done.")

Moreover, I feel this is not the most simple or efficient way to do the bounding box.


texts:
My script works fine, but I am not able to export .pdf file from a report generated. This is my code:

import arcpy

inputFC = arcpy.GetParameterAsText(0)
outputFC = arcpy.GetParameterAsText(1)
matricula = arcpy.GetParameterAsText(2)

query = """ "MATRICULA" LIKE '%s'""" %matricula

arcpy.env.workspace = "C:/Prueba/arcpy"
arcpy.MakeFeatureLayer_management(inputFC,"Tramos_LAT_layer")
arcpy.SelectLayerByAttribute_management("Tramos_LAT_layer","NEW_SELECTION",query)
arcpy.CopyFeatures_management("Tramos_LAT_layer",outputFC)

# Following I try to export report from outputFC

# RLF
plantilla = "PLANTILLA_INFORMES.rlf"

# PDF
pdf = "zz_salida_report.pdf"

# Exportamos Report
arcpy.mapping.ExportReport(outputFC,plantilla,pdf)

I also tested another option, I created a layer file, but it didn´t work either:

lyr = arcpy.mapping.Layer(outputFC)

arcpy.mapping.ExportReport(lyr,plantilla,pdf)

I have checked all post about this issue, although I don't understand why my script doesn't work.


texts: Clipping a SpatialPoygonsDataFrame to another SpatialPolygonsDataFrame

texts:
How do I restrict a SpatialPolygonsDataFrame so that none of the polygons lie outside of an "underlying" SpatialPolygonsDataFrame?

The application I have in mind involves Voronoi Polygons a la Carson Farmer -- the result of such a decomposition is (typically) by default rectangular, but often it's more visually appealing to conform the result to an underlying geometry (city/county/country border, etc.). Another common application is to delete Alaska and Hawaii from US-wide geometries when only the lower 48 are of interest.

This is similar to this question, but more complicated, since both input geometries consist of multiple polygons -- a simple gIntersection will result in an unnecessary dissolution of the main geometry of interest.

It's also similar to this question but I couldn't follow exactly what was going on there because there was no reproducible example. I think the following exercise likely subsumes this question.

Here's a simple example to work with:

library(sp)
library(rgeos)
library(maptools)

points<-list("A" = c(0, 0),"B" = c(1/2, 0),"C" = c(1, 0),
             "D" = c(0, 1/2), "E" = c(1/2, 1/2), "F" = c(1, 1/2),
             "G" = c(0, 1), "H" = c(1/2, 1), "I" = c(1, 1),
             "J" = c(0, -1/2), "K" = c(3/2, -1/2), "L" = c(7/24, 3/8),
             "M" = c(5/8, 3/8), "N" = c(3/8, 5/8),
             "O" = c(5/8, 17/24), "P" = c(3/2, 1))

pts_poly <- function(x)Reduce(rbind, lapply(x, function(y) unlist(points[y])))

poly1_list<-
  list("poly11" = list(ID = "A", coords = pts_poly(c("A", "B", "E", "D" ))),
       "poly12" = list(ID = "B", coords = pts_poly(c("B", "C", "F", "E"))),
       "poly13" = list(ID = "C", coords = pts_poly(c("E", "F", "I", "H"))),
       "poly14" = list(ID = "D", coords = pts_poly(c("D", "E", "H", "G"))))

poly2_list<-
  list("poly11" = list(ID = "a", coords = pts_poly(c("J", "M", "L"))),
       "poly12" = list(ID = "b", coords = pts_poly(c("J", "K", "M"))),
       "poly13" = list(ID = "c", coords = pts_poly(c("M", "K", "P"))),
       "poly14" = list(ID = "d", coords = pts_poly(c("M", "P", "O"))),
       "poly15" = list(ID = "e", coords = pts_poly(c("L", "M", "O", "N"))))

poly1<-SpatialPolygons(lapply(
  poly1_list, function(pl){
    with(pl, Polygons(list(Polygon(coords)), ID = ID))}))
df1 <- data.frame(var1 = rnorm(length(poly1)), var2 = rnorm(length(poly1)),
                  row.names = sapply(poly1@polygons, function(x) x@ID))
polydf1 <- SpatialPolygonsDataFrame(poly1, df1)

poly2<-SpatialPolygons(lapply(
  poly2_list, function(pl){
    with(pl, Polygons(list(Polygon(coords)), ID = ID))}))
df2 <- data.frame(var1 = rnorm(length(poly2)), var2 = rnorm(length(poly2)),
                  row.names = sapply(poly2@polygons, function(x) x@ID))
polydf2 <- SpatialPolygonsDataFrame(poly2, df2)

Lots of code. Here's what we have visually:

(code for plots:)

par(mfrow = c(1, 3))
par(mar=c(0,0,0,0))
par(oma=c(0,0,0,0))
plot(polydf1, xlim = c(0, 1.1), ylim = c(0, 1.1))
title("Underlying Polygons", line = -15, cex.main=1.8)
text(Reduce(rbind, points[1:9]) + .03,
     labels = names(points)[1:9], col = "red")
text(matrix(c(.25, .25, .75, .25,
              .75, .75, .25, .75),
            ncol = 2, byrow = TRUE), cex = 2,
     labels=sapply(polydf1@polygons, slot, name = "ID"))

plot(polydf2, xlim = c(0, 1.6), ylim = c(-.5, 1.1))
title("Polygons to Be Clipped", line = -15, cex.main=1.8)
text(Reduce(rbind, points[-(1:9)]) + .03,
     labels = names(points)[-(1:9)], col = "red")
text(coordinates(polydf2), cex = 2,
     labels=sapply(polydf2@polygons, slot, name = "ID"))

plot(polydf2)
plot(polydf1, add = TRUE)
title("Overlay of Polygons", line = -15, cex.main=1.8)

What I want is to eliminate the parts of triangles E, F, G and H which lie outside of the the the first set of polygons, A, B, C, and D.

A simple intersection is incorrect:

plot(gIntersection(polydf1, polydf2, byid = TRUE),
     main = "Simple Intersection: Overkill")

Here's what my desired output looks like instead (basically eliminating the borders of A, B, C, and D):


texts: How to create GeoTIFF from image and UTM Coordinates

texts:
I am a beginner in GeoTIFF and georeferencing.

I have an image in .emf file format.
 I also have this UTM information: Easting,Northing,Zone,Hemisphere

Example: 400000,4379520,24,N

(I have also the size of the terrain in meters)

How can I use GDAL to create a GeoTIFF from this information?

(My GeoTIFF needs to be in WGS84, EPSG:4326)
(I need to use command line tools, no GUI)

I think I need to convert .emf to a file format which can be supported by GDAL such as .bmp or "standard" .tif. But After ?
Maybe with gdal_translate and/or gdalwarp, but how?


texts: Cannot Edit Features of PostGIS Layers in OpenGeo GeoExplorer

texts:
I admit I  originally posted the question here:
http://getsatisfaction.com/opengeo/topics/cannot_edit_features_of_postgis_layers

But due to no response and not finding a solution, I thought I'd try to ask here.

I have imported some layers from my PostGIS DB into the OpenGeo Suite Geoserver. For the most part they work fine, yet I am unable to edit their features. Every time I try to do this I get the following in FireBug:

feature is undefined
layer.drawFeature(feature,feature.styl...stFeature=layer.selectedFeatures[0];}
OpenLayers.js (line 815)

However, if I import a layer that came bundled with the suite such as World - Borders, everything works fine and I am able to do feature editing as expected. I have also seen this thread:

http://getsatisfaction.com/opengeo/to...

But still I am not sure how to go about it. It is indeed true that we had PostgreSQL and PostGIS installed already, and when I installed the OpenGeo suite all the data dumped itself onto my already existing PostGIS (at least I think so), as everything seems to be running on the same port 5432. What can I do to have my PostGIS layers allow feature editing? This was the main reason we decided to go with GeoExplorer ;).

Geoserver error logs report:

INFO [geoserver.wfs] -
Request: getFeature
handle = null
service = WFS
version = 1.1.0
baseUrl = http://10.0.104.14:8081/geoserver/
providedVersion = null
extendedProperties = {}
query = [net.opengis.wfs.impl.QueryTypeImpl@431e34c7 (group: null, propertyName: null, function: null, filter: [ tib_villages.fid-5a164ff3_134c8d54c27_-6fb2 ], sortBy: null, featureVersion: null, handle: null, srsName: EPSG:900913, typeName: [{http://opengeo.org/#tibet}tib_villages])]
maxFeatures = 20
outputFormat = text/xml; subtype=gml/3.1.1
resultType = results
traverseXlinkDepth = null
traverseXlinkExpiry = null
formatOptions = {}
metadata = null

Thanks for your time and I am looking forward to hearing any feedback,

elshae


texts: Poor Temporal Raster Performance over REST

texts:
When trying to consume temporal rasters (float values) served from ArcGIS Server over a REST endpoint, the results served up are "wonky".  One or two time steps will deliver data as expected but in general, it crashes and returns missing data or inverted values, etc.

If I consume the data in ArcMap directly from SDE, it behaves correctly.  If I consume the data in ArcMap via REST, it is a crap-shoot on what I'll get back, sometimes missing data.

The data is symbolized using "Classified" with the data broken into 7 categories with zero values excluded.  I am currently in the process of converting to integer-based rasters and also reclassifying to see if the issue is corrected but the end-goal is to put this into a web-viewer that I can also query against to get depth values at the x/y click.

The server environment is as follows:

Two Separate Servers...

Windows 2008 R2 - ArcGIS Server 10.1 SP1
Windows 2008 R2 - ArcSDE 10.1 SP1 with Microsoft SQL Server 2012 R2

PS: If I consume the raster using stretched symbology, it returns the data as expected although I can't control the color returned by the server after posting to AGS.  I've tried using FGDB's and registered folders but achieved the exact same issues with that setup.

Thoughts?


texts: How do I calculate ECEF X,Y,Z velocities for GPS if I have no satellites

texts: I'm attempting to simulate a GPS message where know variables are starting location, ending location, speed (m/s), and a bearing (0 - 359 degrees).  How can I calculate the ECEF X, Y, Z velocities with this information when I am not tracking any satellites.


texts: Using satellite images for manual delineation - is pre-processing required?

texts:
I'm new to using satellite images. I'm using Landsat images from USGS (Level 1 data), so terrain-corrected. I want to manually delineate different forest classes using false color composites.

Do these satellite images need to be pre-processed prior to doing this?

I won't be doing any calculations and I don't need pixel values to show reflectance or radiance.

I'm simply wondering whether the images not being atmospherically corrected for will impact their visual look?

I don't want to do all this work if it isn't actually necessary for my purposes, so I'm wondering if someone can help me with this.


texts:
I have been using JavaScript for about a year and a half now, but I still have a fundamental misunderstanding of Promises and can never get them to work.

Can someone tell me what is wrong with this code and why it doesn't wait for the promise to fulfill?

let sitesGeo = [];
let p = new Promise((resolve, reject) => {
     sitesGeo = getSitesGeo(currentBatch, sitesUrl, sitesGeo);
     resolve(sitesGeo);
 }).then((sitesGeo) => {
      console.log("sites geometry ", sitesGeo);
  });

I keep getting undefined for "sites geometry" and I know it's because it is being logged before the promise gets resolved. I know that the function getSitesGeo works because I have another log statement inside of it and it returns the result.

What am I getting wrong, and how to I get console.log statement to wait for the promise to resolve before executing?

Here is the code so maybe this will help to understand why it is not working;

function getSitesGeo(currentBatch, sitesUrl) {
    console.log("getting sites geometry");
    var sitesGeoArray = [];
    var queryTask = new QueryTask(sitesUrl);
    var batchQuery = new Query();
    batchQuery.where = "Batch = " + currentBatch;
    batchQuery.outFields = ["*"];
    batchQuery.returnGeometry = true;
    queryTask.execute(batchQuery).then(function(result) {
        console.log("sites geometry ", result);
        let features = result.features;
        features.forEach(function(feature) {
            var geometry = feature.geometry;
            if(isNaN(parseFloat(geometry.x)) === false) {
                sitesGeoArray.push(geometry);
            }
        })
        console.log("sites array ", sitesGeoArray);
        console.log("geometry type = ", typeof(sitesGeoArray));
        return sitesGeoArray;
    });
}

function getBatchNumnber() {
    var sitesUrl =
    "https://ocean.floridamarine.org/arcgis/rest/services/.../MapServer/0";
    var queryTask = new QueryTask(sitesUrl);
    var batchQuery = new Query();
    var statDef = new StatisticDefinition();
    statDef.statisticType = "max";
    statDef.onStatisticField = "Batch";
    statDef.outStatisticsFieldName = "CurrentBatch";
    batchQuery.where = "1=1";
    batchQuery.outFields = ["*"];
    batchQuery.outStatistics = [statDef];
    batchQuery.returnGeometry = true;
    queryTask.execute(batchQuery).then((result) => {
        console.log("batch query result ", result);
        var currentBatch = result.features[0].attributes.MAX_Batch;
        window.currentBatch = currentBatch;
        //var sitesGeoArray = [];
        async function createSitesGeoArray(currentBatch, sitesUrl) {
            return new Promise((resolve, reject) => {
                const sitesGeo = await getSitesGeo(currentBatch, sitesUrl);
                resolve (sitesGeo);
            });
        }
        createSitesGeoArray(currentBatch, sitesUrl)
        .then(result => {
            console.log("sites geometry = ", result);
        })

I re-wrote the async function to :

async function createSitesGeoArray(currentBatch, sitesUrl) {
            let promise = new Promise((resolve, reject) => {
                const sitesGeo = getSitesGeo(currentBatch, sitesUrl);
                resolve (sitesGeo);
            });
            let result = await promise;
            //return result;
        }
        createSitesGeoArray(currentBatch, sitesUrl)
        .then(result => {
            console.log("sites geometry = ", result);
        })

I don't get any errors, result is now undefined (not waiting to resolve)


texts: How to make the feature clipping by bbox fast in postgis

texts:
I am using postgis to serve 10+ layers, and the data will be returned by a given bbox and a zoom level.

This is the sql I used for query features:

public List<Feature> getFeature(double[] bbox, int zoom){
    List<Feature> features=new ArrayList();
    for(String layer:layers){
        String sql = "SELECT name, ST_Intersection(layer.geom,bbox.geom) as dest_geom FROM layer JOIN (select ST_MakeEnvelope(?, ?, ?, ?, 4326) as geom) bbox ON st_intersects(table_name.geom, bbox.geom) where layer.zoom <= ?";
        ResultSet rs=conn.executeQuery(sql); //

        while(rs.next()){
            Feature f=new Feature();
            f.name = rs.getString("name");
            f.geometry = f.getObject("dest_geom")
        ....
            features.add(f);
        }

    }
    return features;
}

Now it will take almost 200ms for a bbox which conver a 512x512 tile in zoom 16, however it will take almost 10s for a bbox which conver a 512x512 tile in zoom 10.

Though the later have a wide range which means more features, however in a small zoom level the visible features is less too, and that is what the where layer.zoom<=? condition does.

And I have created index for each layer by following this docs:

CREATE INDEX layer_gix ON layer USING GIST (geom);
....

Is there anything I can do to speedup the query?

Update( add explain analyze):

select no-geometry fields:


texts: Getting vertices of geometry that type is not known?

texts:
I am using a WFS service that serves Polygon and Multipolygon geometries. But I do not know ehich ones. So I am converting the wfs response to feature.
this.http.get(url, { responseType: "text" })
  .subscribe((response) => {

    const feature = new GML2().readFeature(response);
    /// get vertices ???
  })

How can I get vertices (corner coordinates) of feature geoemtry? We do not know if that is Polygon or Multipolygon.


texts: Can I have extra row in datatable not used for mapping

texts:
Is there a way that I can disregard 1 row in the table when generating a map and its Jenks values?

We have the need to create separate tables and charts from the data we are mapping.  There needs to be an extra row of data in the data-table for the tables/charts.  But this won’t be included in the map or the Jenks values we’re using for a custom legend.


texts: Fill random data attributes to shapefile using FME?

texts:
there are 150 shape layers and I want to fill the random data in attribute table having different types of fields as Character, integer and float, etc. Each layer have 5 to 15 fields.

Is this is possible to fill random data in all layers using fme?


texts: I want to get the osm of a specific city of USA. I can only get the whole state now. But how can I get the osm data for a specific city.


texts: Desktop tool to batch convert Shapefils to CSV with polygon outline coordinates

texts:
Possible Duplicate:
  How to convert .shp to .csv including attributes AND geometry?

Is there a desktop tool that would allow me to take dozens of large (300MB) shapefiles and convert them into CSV files?  My goal is to extract just a few chunks of data from them (ID, name, lat/lon list of polygon outline) for loading into a database.


texts: Conditional statement in gdal_calc with result as float, not boolean

texts:
Many examples on gis.se that refer to gdal_calc give a boolean example.
example 1
example 2

If I need to calculate a conditional statement on a float:

If Raster_Value >= 1024, then 1024; else Raster_Value.

How do I perform that operation in the gdal raster calculator command line

EDIT
I tried radouxju answer and it worked, except for the cell values that didn't.
here are 3 kinds of results I got

Here is the tool gui and the parameters I've entered


texts: Creating Bivariate choropleth map a.k.a. Relationship map in QGIS

texts:
There is one polygon layer with the number of "immigrants" and "emigrants" respectively. I want to visualize those two values together. Diagram with both normalized values is not an option because the place is already occupied with point symbols.

Of course, the output can be achieved by means of the following approach

"immigrants" and "emigrants" are absolute values, therefore the difference between them can be calculated, i.e. delta = "immigrants" - "emigrants"
Decide what is 'negative' and what is 'positive', i.e.

  CASE
     WHEN "delta" > 0 THEN 'positive'
     WHEN "delta" < 0 THEN 'negative'
     ELSE 'remains the same'
END

In this case, I am considering more "immigrants" versus less "emigrants" as 'positive'.

Define the appropriate amount of categories and visualize them with Graduated symbols, 'positive' with green, 'negative' with red, 'remains the same' with grey.

Another solution can probably be obtained with a half-polygon expression where difference() and $geometry() features are applied. And further, each half polygon visualized with its values accordingly.

However, to go above the typical process I was wondering if it is possible to produce a Bivariate choropleth map a.k.a. Relationship Map in QGIS.

My wish is to visualize both values with two graduated colours simultaneously.
Pink stands for "immigrants" and blue stands for "emigrants" accordingly.

How can I adopt this approach in QGIS? Perhaps, someone has already described the workflow of creating a Relationship Map in QGIS?

I know that ESRI's ArcGIS Online has such option, the requisite procedure was described in this article How to Make a Relationship Map in ArcGIS Online.

Moreover, Dr Cynthia Brewer already described the theory of such mapping in the online guide Color Scheme Types and Combinations: Overview.

Searching on the Internet did not give me any noticeable inspiration.

References:

How to Make a Relationship Map in ArcGIS Online
Effectively displaying demographic data on a printed map
Sequential Sequential Color Schemes
How to present three spatial variables in one map? [closed]


texts: Lat Lon Digitize tool not working with MULTIPOINT in QGIS?

texts:
I don't know what is a reason but in QGIS 3.2 "Lat Lon Digitize" tool doesn't work with MULTIPOINT, it works only with POINT layers.

Is it a bug?


texts: Increasing size/height of OpenLayers map/canvas?

texts:
I'm trying to let the Map/Canvas height of OpenLayers 100%, but I'm not able to do it.

I've already tried putting in several places the "height: 100%;" and used the "map.updateSize()" function, but without success. I already tried to use "height: 100vh;" as well, but it stretches the map instead of increasing the canvas.

I send print screen below illustrating what happens:

I would like the height of the map/canvas to be 100% on the page

That piece of code showed on the image:

        <canvas class="ol-unselectable" width="937" height="480" style="width: 100%; height: 100%; display: block;"></canvas>

        <button class="ol-zoom-in" type="button" title="Zoom in">+</button>
        <button class="ol-zoom-out" type="button" title="Zoom out">−</button>

        <button class="ol-rotate-reset" type="button" title="Reset rotation">⇧</button>

<button type="button" title="Attributions">i</button>

Code where I create the map (Explore.vue):

export default {
  ...
  methods: {
    initMap: function () {
      // eslint-disable-next-line
      this.map = new Map({
        target: 'map',
        layers: [
          layerGroupBaseMap
        ],
        view: new View({
          projection: 'EPSG:4326',
          center: [0, 0],
          zoom: 5
        })
      })
    }
  },
  mounted: function () {
    this.initMap()
  },
  ...
}

Source code: https://github.com/rmmariano/executive_frontend/tree/dev-rodrigo

Does someone have any idea?


texts: Adding annotation feature class to mxd with ArcPy?

texts:
Annotation feature classes in personal geodatabases need to be added to mxds. The annotation is not already in the mxd. I tried Make Feature Layer and add layer through arcpy.mapping and it did not work. Esri states that complex features such as annotation are not supported.

Do I have to open 88 mxds and drag in the annotation feature classes?


texts:
I have a broken package that is preventing QGIS from loading python plugins, updating itself, or even un-installing.

The main symptom of this problem is the following message that appears when QGIS is booted:

QGIS works, but none of the python plugins do and I cannot update, reinstall or uninstall QGIS. The paths of action I have tried so far are:

Reinstall PyQGIS (which I believe to be causing the problem) using Synaptic, which leads to this error:

Updating the packages and reinstalling, using `sudo apt-get update

sudo apt-get upgrade`,

This leads to the following error:

dpkg: error processing python-qgis (--configure):
     dependency problems - leaving unconfigured
    No apport report written because the error message indicates it's a follow-up error from a previous failure.
                                Errors were encountered while processing:
     python-qgis-common
     python-qgis
    E: Sub-process /usr/bin/dpkg returned an error code (1)

Following advice from here on broken packages, forced install sudo apt-get install -f and sudo dpkg --configure -a.

The latter option flags the following error, which tells me (more or less) what the problem is, but not how to solve it:

  dpkg: dependency problems prevent configuration of python-qgis:
   python-qgis depends on python-qgis-common (= 1.9.0+git20121201+80fa7aa~precise1); however:
    Package python-qgis-common is not configured yet.
  dpkg: error processing python-qgis (--configure):
   dependency problems - leaving unconfigured
  Errors were encountered while processing:
   python-qgis-common
   python-qgis

How do I solve this so I can have a fully working version of QGIS again?


texts: Calculation of density metric of 0.75 height quantile based on Airborne LiDAR Scanning point cloud data by lidR package

texts:
This is a piece of code that I wrote:

metrics<-function(z){Q75=quantile(z,probs=c(0.75))

aboveq75=sum(z>Q75)

dh75=aboveq75/length(z)

return(dh75)}

DH75<-grid_metrics(lid,~metrics(Z),res = 30)

I want to calculate the proportion of point cloud quantity above 0.75 height quantile in the total point cloud quantity based on 30-meter grid scale. This metrics can be used to participate in the construction of a forest inventory attributes prediction model using an area-based approach.Because it reflects the distribution of branches and leaves in the upper and middle canopy.But the number of points above the 0.75 quantile is by definition always going be approximately 0.25.Is my definition wrong? or you have a better idea about the Calculation of density metric of 0.75 height quantile.

The result of my definition:

The corresponding value of the grid:


texts: Create Polygon - How do I bring a table in with Lat/long

texts:
I found the below code for creating a polygon. I know what a polygon is, but I'm confused about the code. I have my polygon points defined, but where do I bring my table in to match to the polygon lat/longs? I think it is where I "Set @p = geography", but I'm bringing in a file, not setting a column to a specific lat/long.

IF OBJECT_ID ( 'dbo.SpatialTable', 'U' ) IS NOT NULL
    DROP TABLE dbo.SpatialTable;
GO

CREATE TABLE SpatialTable
    ( id INT IDENTITY (1,1),
    GeogCol1 GEOGRAPHY,
    GeogCol2 AS GeogCol1.STAsText() );
GO

DECLARE @geom GEOMETRY = 'POLYGON(( 55.297622963050465 -4.627166781574488,
                                    55.25851203752759 -3.5724792815744877,
                                    55.034056344339206 -3.5450134612619877,
                                    54.955277168321636 -3.7647400237619877,
                                    54.92372217917785 -4.113555941730738,
                                    54.96946883538248 -4.561248812824488,
                                    55.151935308382306 -4.712310824543238,
                                    55.297622963050465 -4.627166781574488))';

DECLARE @geog GEOGRAPHY = @geom.MakeValid().STUnion(@geom.STStartPoint()).STAsText()

INSERT INTO SpatialTable (GeogCol1) VALUES (@geog)

DECLARE @g GEOGRAPHY
SET @g = (SELECT GEOGRAPHY::STPolyFromText(GeogCol2, 4326) FROM dbo.SpatialTable); -- 4326 is a common SRID

DECLARE @p GEOGRAPHY
SET @p = GEOGRAPHY::STPointFromText('POINT(55.10029008969451 -4.100017953449488)', 4326); -- SRID MUST match polygon

SELECT @p.STIntersects(@g); -- 1 is true, 0 is false

IF OBJECT_ID ( 'dbo.SpatialTable', 'U' ) IS NOT NULL
    DROP TABLE dbo.SpatialTable;
GO


texts: How to display selected features in a map to exporte as PDF

texts:
I've already managed to select feature by ID in a layer and zoom to it using pyqgis. The canvas considered this changes but the map that I exported as PDF do not consider this changes, any idea to solve this problem. This is my code:

xpr = QgsExpression("\"ID_PARCELL\"=1440")
it = vlayer.getFeatures(QgsFeatureRequest(expr))
ids = [i.id() for i in it]
vlayer.setSelectedFeatures(ids)
box = vlayer.boundingBoxOfSelected()
iface.mapCanvas().setExtent(box)
canvas = qgis.utils.iface.mapCanvas() ###
canvas.zoomToSelected(vlayer)
iface.mapCanvas().refresh()


texts:
I have a spatial Grid dataframe, Elev, that I want to crop by a spatial polygon object (state.sp).

I assigned the CRS to both objects with spTransform and the CRS defined as:

CRS.new<-CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")  #EPSG:102003

Both objects seem to be in conformity CRS, because when I type

 > Elev@proj4string
    CRS arguments:
     +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83
    +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

  > state.sp@proj4string
    CRS arguments:
     +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0+datum=NAD83
    +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

Nevertheless, when I do

Felev <- Elev[state.sp, ]

It returns the following error:

Error: identicalCRS(x, y) is not TRUE

Does anyone know why this might be happening?


texts: Count pixels per class instance of a raster classification

texts:
I have a fairly easy problem but I can not find any algorithm for it.
I have this kmeans classification of a windspeed pattern from a model output:

What I want now is to count the pixels in each patch of the classification. Basically meaning I want to count the coherent pixels of each instance. It would be furthermore awesome to have an algorithm to describe the shape of the instances.
What I have found so far is a CNN which does an instance segmentation.
OpenCV instance segmentation

Now I wonder if there is something out there which does the job easier than a NN.
Any algorithm/language is welcome!

Here is the dataset which is plotted:
Dropbox Link to .npy file


texts: Exporting labels from QGIS shapefiles to AutoCAD with labels

texts:
I have shapefiles in QGIS. I want to import labels also in AutoCAD drawing.
I tried to export DXF and open with autocad but labels not visible.
How to import labels in AutoCAD?


texts: Making sense of coordinate systems in retrieved data from a CSW

texts:
I am invoking a CSW server (Catalog Service for the Web) and I am trying to draw the bounding box of each record in the front end of my application.

One of the problems is that the coordinate system of the bounding box is not always in EPSG. What I recieve is something like this:

<ows:BoundingBox crs="urn:x-ogc:def:crs:EPSG:6.11:4326" dimensions="2">
    <ows:LowerCorner>-1.7 -10.3</ows:LowerCorner>
    <ows:UpperCorner>5.4 5.1</ows:UpperCorner>
</ows:BoundingBox>

So the coordinate system information is condensed in a string like this:

urn:x-ogc:def:crs:EPSG:6.11:4326

But this string does not always make sense. For instance, sometimes I recieve strings like this:

urn:ogc:def:crs:::Lambert Azimuthal Projection

Is there a way to transform these coordinate system identifiers into WKT (well known text)? Doing so one can programmatically assign a coordinate system to these strings.


texts: Portal And ArcSDE Geodatabase: Feature Datasets or Not

texts:
We recently had a health check with an Esri representative, to look at our Portal system.  After looking at our ArcSDE enterprise geodatabase she said it was organized very well with everything put into its appropriate feature dataset.  However, she then proceeded to tell me that she does not place feature classes inside feature datasets in cases such as ours, where we use feature layers in a federated server but DON'T use the data store.  She stated that either method works, but she prefers not to put them in feature datasets in cases like ours because how the published features access the database.

We are moving from a test environment to a production environment, and I will be replicating our data over.  Is it best practices not to put data into feature datasets and have them at the root level, will this cause issues with version edits, will it provide a benefit or ruin or data organization scheme?

-Portal - Federated Not Using the Data Store
-Version 10.6.1
-Feature Layers are published and reference the enterprise geodatabase and are NOT hosted


texts: Calculating Cumulative Field Values in ArcMap

texts: In ArcMap 10.2, I wish to create a field which contains the cumulative values of another field in the same attribute table. For example, if the first three cells of the original field contained the values 1, 2, and 3, I would like my new field to contain the corresponding values of 1, 3, and 6. I have seen some solutions to this issue, but as I am not adept in coding I was unable to follow.


texts: Leaflet - resize marker at same time as map zoom animation

texts:
I have a marker (L.divIcon) object in Leaflet. The html content of the marker is an svg. I'm looking to have the divIcon resize when zooming the map.

I've got the necessary code to resize the icons up and down as needed. The issue I'm having is the correct trigger in time. The 'zoomend' event is too late, as by this time the map has already animated. Resizing icons at this stage would have their changes lag behind the map.

map.on('zoomend', function() {
  zoomIcons(currZoom, map.getZoom());
  currZoom = map.getZoom();
});

The 'zoomanim' and 'zoomstart' events don't inform as to the direction of the zoom (in or out) as far as I can tell since map.getZoom() is not updated yet.

Is there another early hook that could be used or some means of determining whether the user is zooming in or out?

For reference below you will find the zoomIcons and dependent functions

/**
 * multiplies CSS values (width, height etc.) preserving their unit
 * @param  {string} cssVal     T css value like -14px or 75%
 * @param  {number} multiplier The multiplication value
 * @return {string}            The css value multipied.
 */
var cssMultiply = function ( cssVal, multiplier ) {
  var returnval = cssVal.substring(0, cssVal.match(/[A-Za-z]/).index) * multiplier
  return returnval + cssVal.substring(cssVal.match(/[A-Za-z]/).index)
}

/**
 * Changes the size of marker icons based on zoom level change
 * @param  {number} oldZoom The previous Leaflet zoom level
 * @param  {number} newZoom The current Leaflet zoom level
 * @return {null}           this function does not return anything
 */
var zoomIcons = function ( oldZoom, newZoom ){

  var zoomRatio = ( newZoom - oldZoom > 0 ) ? 1.5 : 0.666666667;  // 1.25 and 0.8 also works
  var elements = document.getElementsByClassName( 'leaflet-marker-icon' );

  for (var i=0, max=elements.length; i < max; i++) {

    elements[i].style.width      = cssMultiply( elements[i].style.width, zoomRatio );
    elements[i].style.height     = cssMultiply( elements[i].style.height, zoomRatio );
    elements[i].style.marginLeft = cssMultiply( elements[i].style.marginLeft, zoomRatio );
    elements[i].style.marginTop  = cssMultiply( elements[i].style.marginTop, zoomRatio );

  }
}


texts: Creating multiple polygons from a list of coordinate matrices using sf in R

texts:
I'm trying to create an sf object representing the extents of geometry objects in an existing sf object. I would expect that st_polygon() could create a list of polygon objects directly from a list of coordinate matrices, but I can't do it without introducing a clunky for loop.

Is there a faster and/or more elegant way of doing this?

st_bbox_to_coords <- function(bbox) {

    bbox_num <- as.numeric(bbox)
    coords_mat <- matrix(bbox_num[c(1, 1, 3, 3, 1, 2, 4, 4, 2, 2)], ncol = 2)
    return(coords_mat)
}

create_sf_from_extent <- function(sf_object) {

    bbox_list <- lapply(st_geometry(sf_object), st_bbox)
    bbox_coords_mat_list <- lapply(bbox_list, st_bbox_to_coords)

    # What should work
    # bbox_sfg_list <- lapply(bbox_coords_mat_list, st_polygon)

    # What actually works
    dummy_list <- bbox_coords_mat_list
    for (i in seq_along(bbox_coords_list)) {
        dummy_list[[i]] = bbox_coords_mat_list[i]
    }
    bbox_sfg_list <- lapply(dummy_list, st_polygon)
    bbox_sfc <- st_sfc(bbox_sfg_list)
    bbox_sf <- sf_object
    st_geometry(bbox_sf) <- bbox_sfc

    return(bbox_sf)
}


texts: Draw Line in Layout View with ArcObjects

texts:
I'm trying to draw a line on layout view on ArcMap with C# and ArcObjects. I'm working on this sample. But I have some problems.

1.) This code for  draw line on data view but I need to draw line on layout view(In box frame).I want to see properties of line element when I selected my line Item.

2.) When I click on this line in ArcMap, It become invisible..why is this happening?

This is Code Sample :

protected override void OnClick()

    {
        //drawing a polyline

        Point p1 = new PointClass();
        p1.X = 10; p1.Y = 10;

        IPoint p2 = new PointClass();
        p2.X = 20; p2.Y = 20;

        IPoint p3 = new PointClass();
        p3.X = 35; p3.Y = 15;

        IPoint p4 = new PointClass();
        p4.X = 40; p4.Y = 17;

        IPoint p5 = new PointClass();
        p5.X = 50; p5.Y = 19;

        IPoint p6 = new PointClass();
        p6.X = 60; p6.Y = 18;

        IPolyline polyline = new PolylineClass();
        IPointCollection pointColl = polyline as IPointCollection;
        pointColl.AddPoint(p1);
        pointColl.AddPoint(p2);
        pointColl.AddPoint(p3);
        pointColl.AddPoint(p4);
        pointColl.AddPoint(p5);
        pointColl.AddPoint(p6);

        IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
        IActiveView activeView = mxdoc.ActiveView;
        IScreenDisplay screenDisp = activeView.ScreenDisplay;
        short screenCache = Convert.ToInt16(esriScreenCache.esriNoScreenCache);
        screenDisp.StartDrawing(screenDisp.hDC, screenCache);

        IRgbColor color = new RgbColorClass();
        color.Red = 0; color.Blue = 0; color.Green = 0;

        ISimpleLineSymbol simpleLineSymbol = new SimpleLineSymbolClass();
        simpleLineSymbol.Color = color;
        simpleLineSymbol.Width = 1;

        screenDisp.SetSymbol(simpleLineSymbol as ISymbol);
        screenDisp.DrawPolyline(polyline);

        screenDisp.FinishDrawing();
    }


texts: CSV containing GeoJSON format geometries to DataFrame

texts:
I have a CSV file I open as a DataFrame.
Among the columns is one named geom 1 :

df.loc[:, 'geom 1']
>>>
0       {'type': 'Polygon', 'coordinates': [[[2.826589...
1       {'type': 'Polygon', 'coordinates': [[[2.225689...
2       {'type': 'Polygon', 'coordinates': [[[2.225689...
3       {'type': 'MultiPolygon', 'coordinates': [[[[5....
4       {'type': 'Polygon', 'coordinates': [[[3.933055...
                              ...
2998                                                 None
2999                                                  NaN
3000    {'type': 'Polygon', 'coordinates': [[[5.014937...
3001    {'type': 'Polygon', 'coordinates': [[[4.912995...
3002    {'type': 'Polygon', 'coordinates': [[[4.739631...
Name: geom 1, Length: 3003, dtype: object

#Conditionning the dataframe :
df.replace('None', np.nan, inplace=True)

I would like to convert these strings to geometries in a GeoDataFrame.

After a few resarches I did not find a working solution :

TypeError: Input geometry column must contain valid geometry objects :

df.loc[:, 'geom 1'] = df.loc[:, 'geom 1'].apply(wkt.loads)
>>> ParseException: Unknown type: '{'TYPE':'
WKTReadingError: Could not create geometry because of errors while reading input.

Convert GeoJSON to GeoPandas GeoDataframe :

df.replace('None', np.nan, inplace=True)
geom = [shape(i) for i in df.loc[:, 'geom 1'].dropna()]
>>> AttributeError: 'str' object has no attribute 'get'

What I succeded was to sent the DataFrame to PostGIS, convert the type of the column geom 1 from text to geometry and use the function ST_FromGeoJson to access the geometry.
Since I have many columns geom I want to do this under python in order to merge all the columns geom into a single one named all_geoms.


texts: Problem with showing shapefile in GeoServer

texts:
I have a number of shapefiles from a customer.  I can read them with both ogrinfo and QGIS.

They don't have a .prj file, but I know that the projection is EPSG:3006.

When I add them as layers in GeoServer, I declare the projection to be EPSG:3006 in the box "Declared SRS" by inputing "EPSG:3006". I set "SRS Handling" to be "Force declared".

It makes no difference if I add a .prj file with a projection definition. I can then calculate the bounds successfully.

When I try to preview the layer in OpenLayers, I get an empty map. When I try to load the layer as WFS, I get an empty result.

There are no error messages in the logs as far as I can see. I have tried to change the log level.

If I run ogr2ogr and save them as new shapefiles, they work in GeoServer.
If I export them to a new shapefile from QGIS, they work in GeoServer.

Is there a way I can use OGR or QGIS to show my customer what the problem is?


texts: How to Edit custom help in ArcGIS ModelBuilder?

texts:
I have created a model in ArcGIS ModelBuilder with a certain number of layers specified as parameters. When I launch the model from the tool-box, I would like to edit the information in the right hand panel (and replace "no description"). That way I can give information to the users using the model.

Do you know how I can edit that part?


texts: Why QGIS WFS doesn't update PostGIS line and polygon layers?

texts:
I installed PostgreSQL and PostGIS and created PostGIS layers in GeoServer. But when i want to change geometry in QGIS via WFS it drops an error:

":  Error performing insert: java.lang.String cannot be cast to com.vividsolutions.jts.geom.Geometry"

I've read the type of the geometry isn't OK, eg. 'linestring' and 'multilinestring' How could i solve this problem to update my layers?


texts: Where can I find SRTM DSMs (digital surface model)?

texts:
I'm interested in finding a 30 meter DSM (digital surface model that includes canopy) from SRTM with the hope of creating a 30 meter canopy height model. I know that it can be done because the following studies have:

https://www.researchgate.net/profile/Le_Pierce/publication/223834023_Vegetation_height_estimation_from_Shuttle_Radar_Topography_Mission_and_National_Elevation_Datasets/links/02e7e536253c3d0b74000000.pdf

and
http://www2.fiu.edu/~serp1/projects/jpl/jpl.pdf

However, I can't for the life of me find DSMs, only DEMs.  Does anyone have any suggestions for where to look online?  Or do I need to contact JPL or some other agency directly?  I'm also entertaining the idea of ASTER, but am finding those data to be equally challenging to locate.

Where can I find SRTM DSMs? - ResearchGate. Available from: https://www.researchgate.net/post/Where_can_I_find_SRTM_DSMs?tpr_view=K9laIL6pyPlqgwnezqKIxKlevRuyrUeBAiXV_3 [accessed Jan 14, 2016].

I have cross posted because I figure that not everyone uses ResearchGate - https://www.researchgate.net/post/Where_can_I_find_SRTM_DSMs?tpr_view=K9laIL6pyPlqgwnezqKIxKlevRuyrUeBAiXV_3.


texts: QField does not get the GPS location on Android

texts:
I am using QField on three Samsung Galaxy S5 tab tablets for a field work. I created a QGIS project (with project and layers in the same CRS). I downloaded QField on the tablets and created a QField project from QGIS. I imported the project into the tablet.

When I opened the project in QField. I first could not see my location (no pink pin). Then, I upgraded the Android software on the 3 tablets. On 2 of them, it solved the problem but on another one, it did not solve it. On this one, I can see the pink pin wrongly located on the top left corner of the screen but I cannot zoom on location and it is always on the top left corner. I gave the authorization to access my location in all the cases.

What could be the reason for such a problem?

The tablets are the same and I do not understand.


texts: Clip raster without changing cell alignment and avoid resampling on ArcGIS?

texts:
I trying to clip a raster (either by Clip or Extract by Mask) with a polygon.
The original input raster is not aligned to the north but it is aligned to the satellite track.

I have two question:

How do I clip the input raster using a polygon without changing the cell alignment?
How do I clip the input raster using a polygon without resampling the cells.


texts: How to generate multiple maps for each feature using same scale and displaying only the one feature?

texts:
I am using QGIS 2.18.

I've got a layer consisting of 40+ features (lines of movement per day). And I'd like to automatically create maps that always show exactly the same area. But the difference is that I want to show only one feature per map (one day). Is that possible?

I already tried the atlas functionality, but it doesn't provide displaying only one feature and leaving out the remaining, right?


texts: I have been searching on internet for a quite a while and all the sites I found only gave on-site download option (like USGS EarthExplorer), no APIs. Is there any free or paid rest API for accessing 1 arc-second SRTM DEM in .hgt format? And if there are no APIs available what other options do I have?


texts: Getting runtime error from AddFields using ArcPy

texts:
I'm add fields via arcpy code as follows;
import arcpy
arcpy.env.workspace = "C:\Users\sheartz\Documents\NBspray_1952.shp"
arcpy.AddField_management(
    'NBspray_1952',
    [['area_ac', 'DOUBLE', None, 20, None, ''],
     ['area_ha', 'DOUBLE', None, 20,None, ''],
     ['block_1', 'TEXT', None, 12,None, ''],
...more field entry code...
     ['citations', 'TEXT', None, 8,None, '']])

Error was stated as;
Runtime error
Traceback (most recent call last):
  File "<string>", line 46, in <module>
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 3443, in AddField
    raise e
RuntimeError: Object: Error in executing tool

Not seen this error before, can't find what it means.


texts: How to cast box3d_extent to geometry (Rare problem with Geoserver and Postgis)

texts:
I use Geoserver 2.1.3 and "POSTGIS="1.5.3" GEOS="3.2.2-CAPI-1.6.2" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.6" USE_STATS" on PostgreSQL 9.1

From time to time (well - today I can't get around it anymore) I get the following error when I try to "calc native bounds from data" while creating a new layer in Geoserver.

    java.lang.RuntimeException: java.io.IOException: Error occured calculating bounds at
org.geotools.jdbc.JDBCFeatureSource.getBoundsInternal(JDBCFeatureSource.java:485) at
org.geotools.jdbc.JDBCFeatureStore.getBoundsInternal(JDBCFeatureStore.java:179) at
org.geotools.data.store.ContentFeatureSource.getBounds(ContentFeatureSource.java:370) at
    ...
    Caused by: java.io.IOException: Error occured calculating bounds at
    org.geotools.jdbc.JDBCDataStore.getBounds(JDBCDataStore.java:1137) at
    org.geotools.jdbc.JDBCFeatureSource.getBoundsInternal(JDBCFeatureSource.java:478) ...
    94 more
    Caused by: org.postgresql.util.PSQLException: ERROR: function st_envelope(box3d_extent) does not exist
    No function matches the given name and argument types. You might need to add explicit type casts. Position: 30

The geometry column contains points, the SRID is defined in the *geometry_columns* tables. And it worked before... Anyways.. I will upgrade to PostGIS2 in a few weeks.. For now I need help on how to transform a PostGIS box3d_extent to a geometry.

With that information I want to dirty-quick-fix this with the creation of a local hack-function st_envelope(box3d_extent) which would then call st_envelope(geometry).

Thanks in advance for any answers helping with the cast- and/or base problem ;-)

Steve


texts: Suppress verbose output when something goes wrong with ogr2ogr

texts:
Whenever ogr2ogr is not able to read a file, it seems to produce pages of output:

Unable to open datasource `file.geojson' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> NAS
  -> GML
  -> GPX
  -> LIBKML
...

Is there a way to suppress this incredibly annoying output? My workaround is:

ogr2ogr ... 2>&1 | grep -v '^  ->'


texts: Why my calculated NQPD value is different than sDNA calucated?

texts:
I want to calculate the NQPD value (Network Quantity Penalized of Distance) by my hand, and then compare it with snda calculated. But I found that I can't the value same as what sdna calculated.

My case is as the following:

Assume all the links length of this network is 500 meters. While calcualte the NQPD value for the link A with raidus 600 meters and angle metric, totally 7 links should be considered(please refer the bold links in the picture).

In my understanding, the NQPD Angle R600 =1/90+1/90+1/90+1/90=0.0444.
But the sDNA result is 0.03148. Please refer the follwing sDNA result picture.

How does the sDNA get that result and why?


texts: How to use QGIS advancedpythonfieldcalculator from the console?

texts:
I'm trying to use the field calculator in a script like so:

processing.runalg('qgis:advancedpythonfieldcalculator', output_layer, 'Latitude',1,12,5,'','$y',output_layer2)

Not sure exactly how to put the arguments

ALGORITHM: Advanced Python field calculator
INPUT_LAYER <ParameterVector>
FIELD_NAME <ParameterString>
FIELD_TYPE <ParameterSelection>
FIELD_LENGTH <ParameterNumber>
FIELD_PRECISION <ParameterNumber>
GLOBAL <ParameterString>
FORMULA <ParameterString>
OUTPUT_LAYER <OutputVector>

Particularly, the GLOBAL and FORMULA args.  Docs don't show much on this..
I'm getting a return of an empty file (i.e. no points), but the attribute table has the new field name in it. Any ideas?


texts: Openlayers 3 multiple symbols in the same point

texts:
I have a point Vector layer. There are some cases when for one geographic location there are multiple points which should be designed with different symbols. E.g. in a certain building there are pharmacy, grocery store, book shop etc. and each of them has its own symbol.

What is the best practice in this case?

So far I have played with ol.style.Icon objects anchor property, as I have created a field in my db table where I manually set the anchor value of each of them, and use that in my style function. But it is not too efficient...

Does anyone have a good solution for this?


texts: Does map export using GarminCustomMap plugin for QGIS work on Mac OS?

texts:
When I try to export a map using GarminCustomMap (1.4) plugin QGIS 2.01 on Mac OS 10.9 quits with no messages.

I suspect this might be a Maverick related issue.

Is anyone successfully using the plugin with 10.9?

I have reported the crashes on the plugin issues page.


texts: Export theme and its associated layers in xml file in PyQGIS

texts:
I was trying to export all themes and its corresponding layers information through PyQGIS so that  xml can be read later .There is method in Class: QgsMapThemeCollection, writeXml(self, doc: QDomDocument) but cannot able to write into xml file.
Below code has been tried but cannot able to write.
doc = QDomDocument()
qgsl = QgsProject.instance()
mTC = qgsl.mapThemesCollection()
mTs = mTC.mapThemes()
mTc.writeXml(doc)


texts:
I'm trying to create a ol.interaction.Modify() to adjust a bounding box's edges and position. In OpenLayers 2 I used TransformFeature Control and There were options to configure the transformation. How could the feature be modified in ol3 without changing the shape?
I also want to add vertices as in ol2.
If there's no other way then implementing a new interaction, could you provide some resources on how to do that?


texts: Openlayers 4 store drawn feature as TopoJSON

texts:
How can I get drawn features in Openlayers 4 in TopoJSON format. I can get the features in GeoJSON as shown in below snippet

.................................
var draw = new ol.interaction.Draw({
  source: source,
  type: 'Polygon'
});
map.addInteraction(draw);
.................................
function getDrawnFeatures(){
  var format = new ol.format.GeoJSON();
  /* var format = new ol.format.TopoJSON(); */
  /* below line gives me "Uncaught TypeError: format.writeFeatures is not a function" if i use var format = new ol.format.TopoJSON(); */
  return format.writeFeatures(source.getFeatures());
}

I cannot find find writeFeatures function in ol.format.TopoJSON from the docs.


texts: Delete line vertex from polyline layer using point feature in arcpy

texts:
I am trying to delete extra vertex from polyline layer using points coordination's of second layer. My idea is to create from points layer array of coordination and then loop through each line and compare each vertex coordination with point array. Whenever is equal delete it and recreate line.

import itertools, arcpy

arcpy.env.overwriteOutput = True

line_FC = arcpy.GetParameterAsText(0)
point_FC = arcpy.GetParameterAsText(1)
output = arcpy.GetParameterAsText(2)

#clone input lines to memory for fast processing
temp = arcpy.CopyFeatures_management(line_FC, "in_memory/temp")

points = arcpy.CopyFeatures_management(point_FC, arcpy.Geometry())

for point in points:
    with arcpy.da.UpdateCursor(temp, ["SHAPE@"]) as uCursor:
        for line in uCursor:
            line = line[0]
            diff = line.difference(point)
            if not line.equals(diff):
                parts = diff.getPart()

                if parts:

                    joined = list(itertools,chain.from_iterable(parts))

                    poly_trimmed = arcpy.Polyline(arcpy.Array(joined))
                    uCursor.updateRow([poly_trimmed])

arcpy.CopyFeatures_management(temp, output)

I got error:

  Traceback (most recent call last):   File
  "C:\Users\Marry\Documents\Programowanie\del_vertex_outside_poly.py",
  line 20, in
      diff = line.difference(point)   File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\arcobjects\arcobjects.py", line
  1752, in difference
      return convertArcObjectToPythonObject(self._arc_object.Difference(*gp_fixargs((other,))))
  ValueError:  Failed to execute (Script).

What i am doing wrong?

I am using ArcGIS Pro 1.3.


texts:
Sorry I'm a newbie with CRS

I have a shapefile (points) with coordinates in CRS EPSG:3854 (SRP) so that I can show the points within the open layers plugin.

Now I need to transform those points in EPSG:4326. The points do now have Long/Lat coordinates.

But i want to have them in E: and N:

e.g.

LAT: 53.376775 LONG: 9.288940

equals

E:3519305.232 N: 5916149.281

How can I transform them if I have many (like 100) points? Is there a fast way to do this?


texts:
I'm using r.sun in GRASS 7.0.3 to determine beam irradiance at a site (~ 56N -4W).
I have a DEM (a local subset of EU-DEM).

I'm using r.sun in mode 2, called from GUI.

Typical logging in command output:

r.sun --overwrite --verbose elevation=EU_DEM+Trees6@PERMANENT aspect_value=120 slope_value=60.0 beam_rad=BeamIrrad_1021_T6_120as_60inc glob_rad=GlobalIrrad_1021_T6_120as_60inc insol_time=InsTime_1021_T6_120as_60inc day=284 step=0.25

Mode 2: integrated daily irradiation for a given day of the year`

I'm having trouble with the aspect value: this is described as degrees running CCW from East (here:r.slope.aspect and grass70/manuals/r.sun.html), with example that 270 = south.

However, I have undertaken multiple runs altering the aspect_value by 30 degree increments from 120 to 300, i.e. 330, 300, 270... 150 degrees from North.

On an unshaded hill-top, the results suggest that the peak beam irradiance is achieved for aspect_value = 180, i.e. 270degN (west). This is clearly wrong, and suspicious that it occurs at an aspect value that corresponds to South, if the aspect were specified in degrees from North.

The results:

aspect value: 300   270 240 210 180 150 120

equivalent bearing: 150 180 210 240 270 300 330

irradiance: 1161    1972    3323    4481    4881    4439    3268

So, is aspect_value in GRASS 7.0.3 really degrees CCW from east?

(as an aside, it's interesting that the text qualifying the 'aspect' tag notes that it is 'usually measure in degrees from north'...)


texts:
I have a local Nominatim instance in which I carry out searches within a bounding box. The problem is that for whatever value of limit, I get a maximum of 50 results. I understand that this is the limit set by Nominatim so that the public server does not get overloaded.

Consider this reproductible example (sent to the Nominatim public server). Only the server would be different on the local instance.

http://nominatim.openstreetmap.org/search?q=[pub]&format=xml&bounded=1&limit=50&viewbox=4.7718134,45.7073666,4.8983774,45.8082628

As I am running my own instance, I am wondering, is there a way to raise the maximum number of results of my local Nominatim instance?

A possible solution would be to send queries iteratively, with the exclude_place_id argument. For instance, the second query would exclude the place_id of the first query and give a list of 50 other pubs, an so on until no further match is found. But this is not very eleguant...


texts: Grass GIS: Store vector line attribute from raster layer

texts:
I am using SRTM dem (dem raster) projected in UTM, for delineation of streams (Canales) in my area.
I used these grass gis extensions:

r.stream.extract --overwrite --verbose \
      elevation=dem  \
      threshold=1000   \
      stream_raster=Canales  \
      stream_vector=CanalesV   \
      direction=Direccion

r.stream.order stream_rast=Canales  \
     direction=Direccion  \
     strahler=strahler  \
     shreve=shreeve  \
     horton=horton   \
     topo=topo   \
     hack=hack

Now I need to see in the vector attribute table, for every single channel CanalesV its order already obtained in horton and other rasters too.

How can it be done in Grass GIS?


texts: Help choosing an appropriate GIS solution for the following question

texts:
I am currently researching for an appropriate GIS solution. My requirements are pretty simple:

Ability to create circles, symbols or diagrams with transparancy tied to their geographical location (For example showing units sold and total profit for each location)
Size of said symbols should be proportional with values within
They should be capable of displaying values and if possible location names
Most of these informations including coordinates are located in Excel Files
Possibility to automate this procedure every time or at least create semiautomated profiles

I found quite a few solutions such as ArcGIS and Maptitude but I am not certain if such software is required and if there isn't a simpler alternative.

I was also looking at qGIS but sadly I am quite new to GIS and such a form of data visualisation, so any guidance would be appriciated.


texts: Raster processing in ArcPy seems to "simplify"/"round" my resulting values

texts:
I'm working off a DEM and doing some raster algebra with it (real simple, add this decimal value, then divide by this value). The DEM in question is super high resolution LiDAR, that's both spatially and in terms of vertical resolution. So we're working with even small changes in height being significant for this analysis (like centimeters, or hundredth of a meter considering the raster values are in meters).

The problem is the original DEM has very smooth value changes in this gently up-sloping terrain. The resulting raster though is jarring changes, almost like it's been reclassed. Just a big flat area of 0.3, then it jumps up to 1.5 (for example, this is not meters cause it's after the math). There's no gradient between these jumps in value.

I can't figure out what would be causing the issue. Any code that touches the DEM is below:

dem_extract = ExtractByMask(dem, extent_shp)
output_raster = (dem_extract - a) / b
output_raster.save(output_path)

b and a are float numbers, ~0-2.0 range (set by user). The output is Imagine img. Because of smaller a/b values, the output can range from like -1 to 100+. However if the large range is the problem (doesn't make sense as it's a double-value raster, sure the stretch value will change, but that's symbology dependent), I've also tried putting a raster con in addition to the extract to mask any very high elevations too.

Another possible issue is that the DEM is a HUGE mosaic raster (hence the extract first to work with less data and speed up the script). But I don't know why that would make Arc "round" the results and how I would fix that.

The problem is definitely in clipping/masking the raster. In just ArcMap I tried the Data Management version of Raster Clip, it came with the same "steps"/"bins" as the ExtractByMask result with breaks in the same locations, however, the pixel values were now just 1, 2, 3..

So the Extract By Mask tool, despite outputting decimal results, is using some trick (like int(value*100) then converting back with float(new_value)/100.0) but in the process is losing precision in a raster that has such high precision to begin with. Tried looking for a way to just force the mask/clip tools to work in floats but not seeing anything in the environment parameters.

My current stopgap solution is to convert the extent shapefile to a raster (with value just being any >0 field (careful with FID cause it starts at 0). Then use Raster Con to extract where the extent raster is > 0.


texts: ValueError when instantiating LineString from a mix of tuples and Point objects

texts:
I'm new to Shapely, and tried instantiating a LineString object from a mix of tuples and Point objects like this:

>>> LineString([Point(1,2),(2,3),Point(3,4),(5,5)])

I get the following traceback:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7/shapely/geometry/linestring.py", line 46, in __init__
    self._set_coords(coordinates)
  File "/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7/shapely/geometry/linestring.py", line 70, in _set_coords
    self._geom, self._ndim = geos_linestring_from_py(coordinates)
  File "/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7/shapely/geometry/linestring.py", line 228, in geos_linestring_from_py
    "Input %s is the wrong shape for a LineString" % str(ob))
ValueError: Input [<shapely.geometry.point.Point object at 0x1057570d0>, (2, 3), <shapely.geometry.point.Point object at 0x105e69690>, (5, 5)] is the wrong shape for a LineString

According to lines 19-22 of the test script, it's supposed to be able to handle this mix of different objects. Any ideas if this 'feature' is actually implemented? Or should I be sticking to homogeneous class types when instantiating a new LineString object?

p.s. Running python 2.7.5, shapely 1.3.0, Mac OS X 10.9.4


texts: Raster Calculator in QGIS returning 'nan' values

texts:
I am attempting to add two rasters using the QGIS raster calculator and every output ends up as 'nan'. There were no problems with the raster calculator before trying to add these rasters, so I am stumped as to why it is not working with these rasters.

The code is:

"ForType@1" + "LandType@1"

Both rasters have the same projection and the pixel size are the same. It is a simple addition but any sort of calculation with these two rasters, whether addition or multiplication etc, ends up as no data.

Does anyone have suggestions as to why the raster calculator is not working? I am using QGIS 2.6.1, and other users have suggested updating the version to 2.10


texts: Downloading map tiles/image from website

texts:
There are two Norwegian websites that have satellite imagery I need, https://kart.finn.no/ or https://www.norgeskart.no .

Is it possible to save a certain area when zoomed out so the scale is 200m, but at the resolution when zoomed in so the scale reads 20m?

I've tried to see if the images are in temporary internet files but modern browsers make them difficult to find.


texts: merge two tiles of Worldclim data with gdal_merge.py - Nodata

texts:
I am trying to merge two tiles of the worldclim data in geotiff format, which I want to subsequently work with in the raster package in R.

When I use gdal's gdal_merge.py script: gdal_merge.py -v -o tmean1_06_16.tif tmean1_06.tif tmean1_16.tif

I get a merged file, but it seems as if the oceans have all been set to 0. see pic

So I thought that the no data value is not set correctly when merging. Therefore I looked up the no data value in one of the original tiles with gdalinfo.

It seems the no data value is set to NoData Value=-3.39999999999999996e+38

So i tried to use that value when merging:
gdal_merge.py -v -n -3.39999999999999996e+38 -o tmean1_06_16.tif tmean1_06.tif tmean1_16.tif

But that results in a file that seems to either have lost all information or have a legend that is shifted so far, that every temperature value is set to the max. Everything is a uniform green now.


texts: How to remove the extra extent of a raster?

texts:
I recently did a rasterize in QGIS with our city map (which I clipped from the whole country map before) to produce the first map. Now I need to do the connectivity analysis in another software so I did a rasterize and produced the second map with the "extra extent."

My question is how to remove the extra extent (the black area)? I have tried the clip seems doesn't work here. Are there any steps I did wrong cause this?


texts: Difference between Local offset and invisibility depth viewsheds in QGIS?

texts:
Is there a difference between the Local offset viewshed as described by (Fisher, 1996) and the Invisibility depth viewshed available for QGIS created by Čučković (2016)?

Čučković, Z. (2016) Advanced viewshed analysis: a Quantum GIS plug-in for the analysis of visual landscapes. JOSS 1(4). http://joss.theoj.org/papers/10.21105/joss.00032

Fisher, Peter F. "Extending the applicability of viewsheds in landscape planning." Photogrammetric Engineering and Remote Sensing 62.11 (1996): 1297-1302. http://info.asprs.org/publications/pers/96journal/november/1996_nov_1297-1302.pdf


texts: Setting custom center and zoom level for Zoom World icon in OpenLayers 2?

texts:
OpenLayers Zoom Pan Control has a Zoom World Icon. When it is clicked the map zooms all the way out.

Is it possible to override this behavior and define a center and zoom level to which the map should be reset?

My Base Map is Bing Map


texts: XY Co-ordinate information for Start and end point of Line in vector layer

texts: How can we get XY coordinate for Start and end point of Line in vector layer


texts: Coordinate points not lining up with a preexisting map on QGIS

texts:
I asked this question once before, but I did not receive any feedback on solving on the issue and I am still confused as to what I am doing wrong.

I found a public domain map of Florida counties on Google. All I want to do is add coordinate points from a cvs file to the map. It should be simple, right? Apparently not. I made a cvs file with latitude/ longitude and x,y coordinates, but the points do not line up with the map. Yes, I DID have the same projection for both layers. The counties map had a projection of EPSG 2237, NAD 83/ Florida West (ftUS), so I changed the cvs file to have that same exact projection. I also made negative values for the longitude coordinates for United States. I tried this with various Florida maps to see if the counties map I found was just wrong, but every map I tried had the same issue. I tried various types of x and y coordinates for the points, and nothing was better. This really should be a simple thing to do, so I really would like to know what I'm doing wrong.


texts: Ditch Extraction using LIDAR DTM (1m resolution)

texts:
I'm trying to extract ditches using the LIDAR DTM using the software eCognition but I am actually having a hard time since the paper where I based my workflow entitled "Ditches’ network extraction and hydrogeomorphological characterization using LiDAR-derived DTM in wetlands" doesn't have detailed rulesets that were input in eCognition.

I just want to ask, who among you here have an idea how to extract ditches in agricultural landscapes or maybe in an urban setting perhaps, the ones alongside roads?

Or if ever you haven't done such work or project, who else might have a bright idea that will enable me to extract ditches.


texts: Using Go To X,Y dialog box with UTM coordinates in ArcMap?

texts:
Someone gave me x,y coordinate numbers in Meters, and I put them in the Go To X,Y tool of ArcMap, but ArcMap displayed them in the wrong place.

I checked the bottom on the right corner and that says they are in meters.

Is it possible that ArcMap 9.3.1 gave the wrong location when I entered the UTM 13 NAD 1983 numbers?


texts: QGIS 3.4.5 Attribute table expression in print composer not working

texts:
I just upgraded from QGIS 2.18 to the latest stable release (3.4.5).

I am trying to filter an attribute table added to a layout by using @atlas_pagename and I can see there is a match in the expression builder but my table stays blank as if it can't find a match.

This might be the same issue as in this thread (the same problem in 3.4.4) :
QGIS 3.4.4 Attribute table expression not working anymore but on the bug tracker this issue has been closed and reportedly it was solved and corrected for the next release.

Can anyone confirm it's still broken in 3.4.5? I don't know how to report/reopen bugs but will happily report it if anyone can tell me how.


texts: I have seen other people mention the issue of text annotation scaling which does not match the scale on the map once you convert the document into a pdf. Does anyone know if there is a way to resolve this in QGIS 3.6.1? I have tried to increase the font to a crazy size but it still looks small on the pdf. If I look up the 'Item' menu, there is no option to change rendering as suggested in previous blogs. Any suggestions to help a novice like me?


texts: Working with Sentinel-2 data and aggrigation in R

texts:
So I was wandering if there is a way to simplify a couple of Sentinel's rasters so I can apply raster analysis on them. Problem is that whenever I try to aggregate the data my RAM, CPU and storage overload (see the following line). And that is by using parallel processing.

r2=aggregate(ras,fact=10,fun=mean, expand=TRUE, na.rm=TRUE)

Now I am wandering, if there is a faster way to transform my data, faster, using another method or function. Now, I have also thought of resample and splitting but I wanted to ask you people if you knew a faster way to aggregate my data.


texts: Testing for matching shapefile pair then processing pair Arcpy

texts:
just looking for some help to solve a logic issue (very inexperienced with python/programming).

The task I have is to check two different folders for a matching pair of shapefiles where the strings would match up to the first "_". The inherit problem with this was that some shapefiles would be missing their counterpart and I only want to process if both exist. The only way I could figure out how to make matching pairs was to rename the shapefile with a split and store in a feature list then test for the match.

The problem with my solution is that sometimes the script will need to be re-run on the shapefiles, so I don't really want to rename the data or have to back it up before each run manually (it is part of a statistical cleaning and regression model for the data, sometimes the data needs to be scaled after the fact or additional sources are added).

How do I do this without renaming the shapefiles?

import arcpy, os, numpy, datetime
from arcpy import env

env.overwriteOutput=True

harvest = arcpy.GetParameterAsText(0)
applied = arcpy.GetParameterAsText(1)
output = arcpy.GetParameterAsText(2)

#Create or clean the working gdb
if arcpy.Exists(os.path.join(output, "cleaning.gdb")):
    arcpy.Delete_management(os.path.join(output, "cleaning.gdb"))
    arcpy.CreateFileGDB_management(output, "cleaning.gdb")
else:
    arcpy.CreateFileGDB_management(output, "cleaning.gdb")

#save current time
now = datetime.datetime.now()

#Create applied feature list and rename shapefiles
env.workspace = applied
fcs2 = arcpy.ListFeatureClasses()
print("Renaming As Applied Data")
for fc in fcs2:
    fn = fc.rsplit('_',3)[0]
    print(fn)
    fc.replace(fc, fn)
    arcpy.Rename_management(fc, fn)
fcs2 = arcpy.ListFeatureClasses()

#Create harvest feature list and rename shapefiles
env.workspace = harvest
fcs = arcpy.ListFeatureClasses()
print("Renaming Data")
for fc in fcs:
    fn = fc.rsplit('_',3)[0]
    print(fn)
    fc.replace(fc, fn)
    arcpy.Rename_management(fc, fn)
fcs = arcpy.ListFeatureClasses()

env.workspace = r"path"

#Test to make sure we can match!
print("Checking for matching datasets!")
count = 0
for fc in fcs:
    if fc in fcs2:
        #if there is both data for a field we will clean!
        #set variables first for file paths
        adata = os.path.join(applied, fc)
        hdata = os.path.join(harvest, fc)
        outputtable = os.path.join(output, fc.rsplit('.',1)[0] + ".xls")
        print("> We have a match!")
        print("> " + fc)
 ##run method etc


texts:
I've been started developing GoogleEarth API on web. I'm using GoogleEarth plugin 7.1.2 version. It is probably a simple question but I'm having trouble with that.

I'm trying to integrate the ArcGIS Image Extension for Server (version 10.1) to a web application based on google plugin and server enterprise, or at least find a proper solution. I know that there is a way to publish ESRI service as KML but I want it to be a last resort cause the image service is outsource

Sorry if my English is not that good.


texts: QGIS: exporting contours to DXF with height values

texts:
I am using this video tutorial to create ground contours for a noise model in iNoise.
https://www.youtube.com/watch?v=JxCxBREEqG0
See tutorial 8 https://dgmrsoftware.com/products/inoise/gallery/
When I import the dxf file it has no height data.  Not sure where I'm going wrong as all the layers have height data when checked in QGIS.
It seems I loose it when exporting as a dxf.
This question may indicate that dxf has no height values but is drawn to height.  Not sure.
Exporting QGIS contour attributes in DXF file?


texts: Convert EPSG:2263 LiDAR data to EPSG:4326 to use with FSX

texts:
I have LiDAR (las) input files in EPSG:2263 (NAD83, NY state plane, feet_US) aquired from https://gis.ny.gov/elevation/ and want to convert those file to EPSG:4326 (WGS84) for use in MS Flight Simulator. I tried gdal / pdal/ liblas, lastools etc. but it seems impossible to do the right conversion. I've seen a similar question (Using proj to convert EPSG:2263 to WGS84) but there are only pairs of x and y coordinates converted. I can't see an option for z (Height) coordinates.

I have tried to create a delimitered text file first with las2txt, and then i used gdaltransform -s_srs EPSG:2263 -t_srs EPSG:4326 <in.txt >out.txt. It produces a text file with lat long coordinates, but the conversion is somehow wrong (see image).

What is the proper way to do the conversion?


texts: Passing Python List into Layer Definition Query

texts:
I am trying to cursor search through items in a field in a feature class. I then append them to a list called 'distList'. I want to pass the list into layers definition query as follows:

"Field1" in distList

or in long hand, it would look like this:

"Field1" in ('ds(c1)', 'ds(b1)', 'ds(c2)', 'ds(g1)')

The list looks like this:

['ds(c1)', 'ds(b1)', 'ds(c2)', 'ds(g1)']

The problem I am having is, the definition query above doesn't like to look 'in' a python list. It doesn't recognize the square brackets. It prefers round brackets, so I thought about using a tuple. The problem I am facing is I don't know how to build a tuple after cursor searching the items from the feature class. here's my code example:

distList = []
for lyr in arcpy.mapping.ListLayers(mxd, "", df):
    if lyr.name == "Disturbance":
           for row in arcpy.SearchCursor(lyr):
                if "ds" in row.Field1:
                     distList.append(row.Field1)
    lyr.definitionQuery = '"Field1"' + "in " + distList

Can anyone suggest a way to get my list into a tuple or maybe just a better way to get my items into a format where they have round brackets instead of square ones?

As a work around, I tried converting the list to string str(distList) and then replacing the brackets. That works fine, but it seems cumbersome and I'm sure there's a better way.

Thanks, Mike


texts: Export data driven pages indexed by a string value

texts:
Below is a portion of a script that works with a ddp enabled map document using an index grid where the name is a string (examples A23, B34, etc) and the sort field (ie the page names) are simple integer values. The first portion makes a list of page numbers based on a query and the second part exports those pages as PDFs.

I'm trying to modify it for a different index grid where the sort field is the same field as the string name field so the range function doesn't work. I can get a list of all the page names that I want to export but I can't figure out how to acutally get them to export. I figure its something along the lines of "for items in names_list:...." but I can't find a list of ddp functions anywhere so I'm really stuck. Anyone know how to make this work or alternatively does anyone know where a list of ddp python functions and their use exists?

for row in arcpy.SearchCursor(grid):

    expression = str(row.getValue(field_name))
    if expression.startswith(letter_block):
        x = row.getValue(field_num)
        range_list.append(x)

    for page_num in range(range_list[0], range_list[len(range_list)-1]+1):
        mxd.dataDrivenPages.currentPageID = page_num
        page_name = mxd.dataDrivenPages.pageRow.getValue(field_name)
        arcpy.mapping.ExportToPDF(mxd, str(page_name) + ".pdf")


texts: Scrolling through selected records in attribute table

texts: I am working with data from GPS tags to determine departure/return times so am I am currently going through the attribute table for each tag and clicking each consecutive record to highlight the point on the map. Is there a way to continuously "scroll" through selected records in the attribute table? I tried with a mouse wheel and arrow keys but they just select one cell of the attribute table and it doesn't highlight the whole record so the point is not highlighted on the map.


texts: SVG to GeoJSON conversion - incorrect render

texts:
I have a vector image in SVG format.(without georeferences)

I am trying to display it using OpenLayers to display properties of different polygons in it.

I tried converting it to GeoJSON but the resulting render appears to be stretched (ratios are not preserved).

Is there a transformation or a map projection that I am missing  or recreating a GeoJSON from scratch the only relevant choice?


texts:
I had been very glad if any of you have time to help me.
I am a student and currently writing a thesis.
In this paper I will create a 3D model of a glacier, where I will show the internal structures of the glacier. (Ice and Snow) I use ArcMap and ArcScene.

My workflow is as follows:
Collected data with GPR (ground penetrating radar)
Analyzed data and created given depth values (how thick the snow is and how thick is the ice). This is done in Reflexw. The output I get is in an ASCII format, columns of snow thickness and the thickness of the ice.
Of this data I create shapefiles in arc map. A Shapefile for snow and the ice. I have about 7 sections with snow and ice thickness.
My goal now is to create one TIN for the snow layer and one for the ice layer. I already have a TIN of the research area with orthophotos draped over (topography). Now I want snow and ice TIN lie below each other, and the TIN model of the topography. Have tried to create negative values of the depth measurements, but to no luck.

My question is: How do I create TIN model for snow and ice, so that the layers lie below the topography? :)

I'm attaching some links so that you may understand better what I mean:

This is what i want to do:

http://resources.arcgis.com/en/help/main/10.1/index.html#//00q900000032000000

and now it looks like this: http://imgur.com/yYscxbp (the blue and red lines is the sections of snow (red) and ice (blue)). The green TIN is made from the snow layer, but as you can see; it does not fit with the topography. Do I have to adjust the z coordinatesystem in some ways?

Highly appreciate all the help I can get.

Regards, Jorn


texts: Is there a way for multiple users to use an SQLite/SpatiaLite db simultaneously?

texts: Is there a way that multiple users can use the SQLite / SpatiaLite database simultaneously?


texts: Make Feature Layer / Select / Select Layer By Attribute "Get Unique Values" grayed out in ArcMap ModelBuilder [10.8.1]

texts:
I'm trying to extract road segments from a feature class through ModelBuilder in ArcMap 10.8.1 based on attribute queries, but the "Get Unique Values" button does not pull up a list of values when I click on it. Of the three tools that I have tried for this (Make Feature Layer, Select, Select Layer By Attribute) they all produce the same non-result.
Below is a screenshot that shows the tools I've used leading up to the selection step, and what happens when I try to "Get Unique Values" from these various tools. I know that sometimes ModelBuilder is fussy about which selection tool is used based on whether the data is a feature class or feature layer.
Can someone tell me what is happening and how to make it work?

Cross-posting from GeoNet to try to reach as many people as possible!

UPDATE
I should mention that if I try to manually type a query into the main Make Feature Layer, Select, or Select Layer By Attribute tool, it won't verify.
Still haven't figured out the initial problem, but I did find a workaround by inserting a Calculator Expression from the Create Variable option and linking it to the tool. Hopefully this is helpful!


texts: Adding and populating Geometry column from an existing table

texts:
I have a PostGIS table with three columns (rid, rast, fileName), now I want to add a geometry column and populate the geometry column using the rast column. A single rast index may have multiple geometries.
Normally to get geometry I use:-
SELECT geom FROM elevationData,LATERAL ST_PixelAsPoints(rast, 1)
How can I add and populate the geometry column? I am getting confused on how to use the WHERE clause in this scenario.


texts: Calculate the diameter of an `sf` object

texts:
The diameter of a set of points in a metric space is defined as the supremum of all pairwise distances between points in the set. So in the case of the earth's surface, the diameter of a (closed) set of points is equal to the diameter of the smallest circle that encloses the whole set.

Given an arbitrary sf object, such as a MULTIPOLYGON, how can I compute the diameter?


texts:
I'm trying to find if there is some solutions for collect data from users via some easy form (not google forms, something more powerful)?

I heard about arcgis collector, but you have to log in (so not everyone can submit a form) and it's mostly for organisation rather than for a large public..

Is it possible to collect this data not by manual writing input from users, but by simple click on a map joined to the form?


texts: Query shortest path from coordinates with pgrouting

texts:
How do I calculate the path, pass coordinates and, if possible, in the return function, have the route coordinate list?
(I'm trying to plot the route with the ipyleaflet library in a jupyter notebook, and I need the return of the route's coordinate list. I'm new to this and I don't know how to do it)

I tried with an example query that I searched for in a book
but I got an empty answer

SELECT r.seq, r.node, w.name As street, w.x1, w.x2, r.cost::numeric(10,4), ( SUM(w.length_m) OVER(ORDER BY r.seq) )::numeric(10,2) As dist_m
FROM pgr_dijkstra('SELECT gid AS id, source, target, cost, reverse_cost FROM ways', (SELECT id FROM ways_vertices_pgr ORDER BY the_geom <-> ST_SetSRID( ST_Point(-3.7482, -38.5486), 4326) limit 1 ) , (SELECT id FROM ways_vertices_pgr
ORDER BY the_geom <-> ST_SetSRID( ST_Point(-3.7532, -38.5170), 4326) limit 1 ), true) AS r
LEFT JOIN ways AS w ON r.edge = w.gid;


texts: Accessing and manipulating fields in Layer file

texts:
I am trying to update several fields in a mosaic dataset, to do this I am creating a mosaic layer and then running a loop over each row and calculating fields.
A NAME field contains a string which I would like to use to calculate several other fields. It is formatted like the following:
Project_2020_1x_YZ_Type_20200101_0900VLAT_102a

I would like to use split and use the output list indices in a field calculation (ie Project, Year and Date). The issue is I don't know to access the NAME field within a mosaic layer.
# create mosaic layer
mdl = arcpy.MakeMosaicLayer_management(md_path, "md_layer")

# calculate fields loop
print("Calculating fields in footprint...")
for lyr in mdl:
    split_name = lyr.NAME.split('_')
    arcpy.CalculateField_management(lyr, "YEAR", ????, "PYTHON")
    print("Year calculated...")
    arcpy.CalculateField_management(lyr, "PROJECT",????, "PYTHON")
    print("Project calculated...")

It falls over at the first step - split_name = lyr.NAME.split('_')-
AttributeError: 'Layer' object has no attribute 'NAME'

I would like to know the best approach to accessing specific fields within the mosaic layer.


texts: Python script to accumulate flow volumes for a pipeline geometric network

texts:
I have a point wells layer and a polyline pipeline layer. The wells are snapped to the end of the pipeline segments. The wells layer has weekly updated flow volumes and I'd like to find a way to automate the transfer of this flow volume from the point layer to a new field in the line segment that the point is attached to.

(I'm aware of spatial joins but this creates a new feature class which won't work because the pipeline layer is part of a geometric network. I just need to transfer the attribute to the original pipeline layer)


texts: Distribute the points randomly within a polygon in QGIS

texts: Say i have 20 points within a polygon. Each point represents a feature. out of which 15 were distributed and 5 points were one above the other in same location. Is there any way to randomly distribute those five points using QGIS?


texts: UTM coordinates and knowing how to get the grid zone letter?

texts:
I am using a framework called GeodeticUTMConverter to convert from decimal style coordinates to UTM. However the framework only provides the Northing, Easing, Grid Zone and Hemisphere, but not the grid zone letter.

If I have the following:

northing = 5067804.5036825798
easting = 551203.31329645123
gridZone = 10
hemisphere = kUTMHemisphereNorthern

The UTM would be 10T E 551203 N 5067804.

What is the most accurate way to calculate what the grid zone letter should be? How does the T get calculated?

EDIT: Added method that I wrote to calculate:

+ (NSString *)UTMlatZoneForCoordinates:(CLLocationCoordinate2D)coordinate {

    NSString *DigraphLetrsE = @"ABCDEFGHJKLMNPQRSTUVWXYZ";
    int latz = 0;//Latitude zone: A-B S of -80, C-W -80 to +72, X 72-84, Y,Z N of 84
    double latd = coordinate.latitude;

    if (latd > -80 && latd < 72) {
        latz = floor((latd + 80)/8)+2;
    }
    if (latd > 72 && latd < 84) {
        latz = 21;
    }
    if (latd > 84){
        latz = 23;
    }

    char letter = [DigraphLetrsE characterAtIndex:latz];
    return [NSString stringWithFormat:@"%c", letter];

}


texts:
I am totally new to Qgis (or GIS in total):
My problem is to get two layers in same CRS to do an intersection and further analysis described in the blog of underdark (Groupstat tutorial).

My is the CORINE-AT clip from here: http://www.umweltbundesamt.at/umweltsituation/umweltinfo/opendata/oed_landbedeckung/?cgiproxy_url=http%3A%2F%2Fhttpapp5.umweltbundesamt.at%2Fdata%2Fdataset.jsf%3Bjsessionid%3D2E72E8BB3EB1C39EEBDC124E516B445C%3Fid%3D1

and the country data from DIVA-GIS: http://gadm.org/data/shp/AUT_adm.zip

As far as I could see the CORINE data is in Austrian Lambert and the country data is in WGS84.

I tried to change the CRS via right click to layer, via vector-> define actual projection, change the WGS84 to Austrian Lambert and then save the shape file...

There were always two outcomes:
Error in CRS:... or

The layer did not lay on each other (country data is flattened).

However the result when doing the intersecting (or clipping, or anything else) is an empty attribute table.


texts:
I have a map in ArcMap that I would like to export as an image, adding also a legend and so on. My issue is that when I go to the Layout View, my data occupy only a small part of the layout view. I would like to zoom in the data, but if I try I do not zoom the data, but the whole layout view.

How can I zoom in in layout view to enlarge the view of my data without zooming the layout view?


texts: CONVERT lat and long field in PostgreSQL using ST_GeomFromText

texts: i have two fields for lat and long in point text record.how can i use ST_GeomFromText function for all of my records to convert it into geometry??


texts: Is there a viewer function/display in QGIS as it is in ArcGIS. Or is there an other possibility to limit the access (just visualisation, no analysis)?