texts: Creating new layer with only subset of polygons of existing layer in GRASS?

texts:
Suppose I have a map with 10,000 polygons with categories 1 to 10,000 in layer 1. I would like to create a second layer with only the features in layer 1 with categories 1 to 5,000. From the online documentation, it seems I can do it with a combination of v.category, v.to.db and v.db.connect. But v.category does not provide a filtering mechanism when creating a new layer.
It seems doable from this link, specifically:

  If a vector object has zero categories in a layer, then it does not appear in that layer. In this fashion some vector objects may appear in some layers but not in others. Taking the example of the fields and paths, only some boundaries, but not all, might have a category value in layer 2.

I googled very hard but did not find a single example so far.


texts: define a datum and meridian in a CRS in QGIS

texts:
I have a shapefile defined by the source as:
"The layer is in geographical (lat-long) coordinates (GCS_Clarke_1866; Datum: D_Clarke_1866; Prime Meridian: 0)."

opening the shapefile in QGIS 1.8 the proj4 format is:

+proj=longlat +ellps=clrk66 +no_defs

Why is there no datum defined and how can I define it? Is it necessary to define the Prime Meridian?


texts: Merge Two layers Keeping only common fields

texts:
I have two feature classes that I need to merge to a single layer. There are ~10 common fields that I need to preserve and each of the two layers has some fields that are not shared.
I know that I could first DeleteField_management() each layer to have layers with identical fields and then do a simple merge but it seems like I should be able to use Merge_management() and a field mapping. I don't see how to include/exclude fields though, only how to transform them.
Inputs: file_1 with fields a,b,c,d. file_2 with fields c,d,e,f.
Desired output: file_3 with fields c,d.


texts: Create geometry column in Geopackage with SQL

texts:
How to create a table with geometry column that is recognized as spatial layer? I use the DB manager from QGIS 3.4.2

I use the following:

CREATE TABLE 'sample_attribute_with_geometry' (
  id INTEGER PRIMARY KEY AUTOINCREMENT,
  geometry MULTIPOLYGON NOT NULL,
  name TEXT NOT NULL,
  description TEXT NOT NULL,
  date INTEGER NOT NULL,
  source TEXT NOT NULL)

The table created is not spatial.

How could i do?


texts:
I've produced a map in QGIS, with various shapefile layers.  I want to share this with a colleague.

If I send him everything, keeping the same file structure, will it open ok on his machine?


texts: "non finite transformation detected" using sp R

texts:
This is a part of a broader question already posted (Non finite transformation detected - lidR spTransform) but I am creating a new question for the sp related part of that.
I'm trying to reproject points using spTransform and get non finite transformation detected error. Below is a reproducible example with which I was helped by @JRR:
library(sp)
X = c(6275235.88, 6275229.48, 6275163.35)
Y = c(2094771.54, 2094815.15, 2094777.31)
WKT = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",39.83333333333334],PARAMETER[\"standard_parallel_2\",38.33333333333334],PARAMETER[\"latitude_of_origin\",37.66666666666666],PARAMETER[\"central_meridian\",-122],PARAMETER[\"false_easting\",6561666.667],PARAMETER[\"false_northing\",1640416.667],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]],VERT_CS[\"NAVD88 height - Geoid12B (ftUS)\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AUTHORITY[\"EPSG\",\"6360\"]]]"
XY = cbind(X,Y)
XY = sp::SpatialPoints(XY, sp::CRS(SRS_string = WKT))
sp::spTransform(XY, sp::CRS("+init=epsg:3857"))
#> Warning in sp::spTransform(XY, sp::CRS("+init=epsg:3857")): 3 projected point(s)
#> not finite
#> non finite transformation detected:
#>            X       Y
#> [1,] 6275236 2094772 Inf Inf
#> [2,] 6275229 2094815 Inf Inf
#> [3,] 6275163 2094777 Inf Inf
#> Error in sp::spTransform(XY, sp::CRS("+init=epsg:3857")): failure in points 1:2:3

I also got the same error when projecting to other CRS like EPSG 6417. I know this must be because I'm asking something mathematically impossible but I don't see what as it looks like all the info are available to perform the reprojection.


texts: What in-depth spatial database systems tutorials exist?

texts:
Is there a good tutorial that explains in depth the internals of a GIS and a spatial database system such as PostGIS with examples and without any background assumption? I'm specifically looking for answers to questions such as:

What is a geometry in GIS?
How is it represented?
Given a lat/long coordinate, what are the operations I need to perform on it to bring it to a state where I can call something like ST_contain PostGIS to check if the polygon contains the lat/long coordinate?

I have some experience with general relational database systems but no background with GIS or the spatial database paradigm.


texts: Set the right crs on sf object to plot coordinate points

texts:
UPDATE
This question has already been answered on stackoverflow.

I'm trying to define the right crs for my sf object. I want to plot points on the following layer (country: the Netherlands):

Simple feature collection with 380 features and 3 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 619232.6
epsg (SRID):    NA
proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs

output:

This layer has the correct projection.

But the POINT layer doesn't have the correct crs project presumably because it has no proj4string?

Simple feature collection with 566 features and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 3.5837 ymin: 50.86487 xmax: 7.120998 ymax: 53.44835
epsg (SRID):    NA
proj4string:    NA

How do I set the same projection as the previous map so I can plot the coordinates points on it?

When I try it in QGIS and select the same crs for both maps they seem to be each on a whole different world.


texts:
I am trying to build a SQL expression using SQL Mode in "Make Feature Layer" tool. I have a Land Use/Land Cover Layer that contains several residential classifications (low density, medium density, high density, rural). I want to build a SQL statement that makes the Feature Layer out of any feature that has the word "RESIDENTIAL" in their label. I can't seem to get the logic right.

I thought it was going to be like ArcMap where they have the "SELECT * FROM WHERE" already implied, but that doesn't seem to be the case.
I tried "RESIDENTIAL" IN LABEL12, like how I would in ArcMap, but that didn't work either.

My statement is as follows:
          SELECT * FROM "Land Use 2012" WHERE LABEL12 LIKE "RESIDENTIAL"


texts: Create multipoint layer from array by expression

texts:
I have two layers. One shapefile point layer, called "point1" and one csv no geometry layer, called "Meta".

"point1" looks like:

id   bd
1    a
2    b
3    c
4    d
5    e
6    f

and so on...

"Meta" looks like:

Name    ids
sample  1, 2, 3, 6

It is a one-line layer and 1, 2, 3, 6 is one single string.
Layer "point1" will always stay the same and "Meta" will change everytime the project is opened.

Now I want to create a new layer (multipoint) from point1 layer which id's are listed in the ids string of Meta layer. Out of curiosity I want to know whether this is possible or not using Geometry by expression. (If you have an answer for virtual layers feel free to share as well.)

So far I have the following expression built up which returns an array like this: [<feature: 0>, <feature: 1>, NULL, <feature: 3>] and so on... Here is my expression:

array_agg(
    get_feature( 'point1','id',
        array_find(
            (
                array_foreach(
                    string_to_array(attribute(get_feature_by_id('Meta',2),'ids')),to_int(@element)
                )
            )
                ,attribute( $currentfeature, 'id' )
        )
    )
)

Replacing array_agg() with geometry() returns a copy of point1 layer keeping only the wanted geometries. But how can I create a multipointlayer out of it?

Can I create a multipointlayer using expressions and arrays using QGIS 3.6 this way?


texts: Transforming coordinates to a topocentric system using PROJ.4?

texts:
I would like to transform coordinates (e.g. WGS84 or Bejing1954) to a topocentric system.

There seem to be several EPSG codes for this (5819, 5820, 5821, 15594, found here or here) but none of them is working with PROJ.4

It just shows the message: no options found in 'init' file

I saw that it is possible to define a custom transformation in PROJ using several parameters, but since I am not familiar with this I have no clue how to do it.

Is this possible at all? And i f yes, could someone help me to compile the correct set of parameters?


texts: How to obtain detailed Storage usage from ArcGIS Online?

texts:
The My Organization tab in ArcGIS Online shows that I've been using lots of credits for Storage.

It's possible that most of this cost is due to a single large dataset that no-one uses, but I can't see any options to sort by size, how frequently the datasets were accessed, etc.

Is it possible to determine which of the many layers in my account are using the most credits?


texts: Can SDEMON command be used on ArcGIS Desktop?

texts:
Working on ArcSDE I need to use SDEMON to monitor or stop users, so that I can do my updates.

Can SDEMON command be used on ArcGIS 10.2 for Desktop?


texts: Rasterizing Shapefile in different projection/extent from sample raster in Python?

texts:
I have a road network shapefile that I'm trying to rasterize and combine with the DEM raster. I created the road network shapefile from a dataset (using fiona) and therefore it has no projection. What I'm trying to do is convert this shapefile to a raster with same projection and extent as the DEM raster. My code runs without nay errors but I get a blank raster with 0 values.

Here's the code I'm using:

def vector_to_raster(vector_path, raster_path, output_path, x_size=None, y_size=None):

'''
This function creates a raster object by burning the values of this
shape file into a raster with the given resolution.
'''

vector = ogr.Open(vector_path)
source_layer = vector.GetLayerByIndex(0)
x_min, x_max, y_min, y_max = source_layer.GetExtent()
print (source_layer.GetExtent())

#x_resolution = int((x_max - x_min) / x_size)
#y_resolution = int((y_max - y_min) / y_size)
#print (x_resolution, y_resolution)

## Reference Raster
raster_ds = gdal.Open(raster_path, gdal.GA_ReadOnly)

# Fetch number of rows and columns
ncol = raster_ds.RasterXSize
nrow = raster_ds.RasterYSize

# Fetch projection and extent
proj = raster_ds.GetProjectionRef()
ext = raster_ds.GetGeoTransform()

# Create the raster dataset
memory_driver = gdal.GetDriverByName('GTiff')
out_raster_ds = memory_driver.Create('Samoa.tif',ncol, nrow , 1, gdal.GDT_Byte)

out_raster_ds.SetProjection(proj)
out_raster_ds.SetGeoTransform(ext)

status = gdal.RasterizeLayer(out_raster_ds, [1], source_layer, options= ['ALL_TOUCHED=TRUE',  # rasterize all pixels touched by polygons
                              'ATTRIBUTE=Shape__Len'] )

#if status != 0:
#    print("I don't think it worked...")
#else:
#    print("Success")

return

Here attribute 'Shape_Len' is what I want as pixel values in the new raster.

As a workaround, I tried to create a raster from the shapefile without giving any reference raster and using shapefile's extent to get the pixel width and length. In this process I got a raster with the values that I wanted but now the two rasters are in different projection and extent. When I try to bring both the rasters to same extent and projection, I again get a blank raster.

Any thoughts on this?

def reporject_raster():
dst_crs = 'EPSG:4326'

with rasterio.open('raster_from_shapefile.tif') as src:
    try:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })
    except CPLE_AppDefined:
        with rasterio.open('reference_raster.tif') as src2:
            transform, width, height = calculate_default_transform(
            src2.crs, dst_crs, src2.width, src2.height, *src2.bounds)
        kwargs = src2.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })

        with rasterio.open('output_reprojected.tif', 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src2.transform,
                    src_crs=src2.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest)


texts: Showing Progress Dialog in ArcGIS Desktop using ArcObjects?

texts:
I have an ArcGIS custom IGPFunction tool.  I want to use the GP progress dialog not the progress bar at the bottom of the ArcGIS 10 window as I don't want people to interact with the map whilst it's processing.

In the Execute method you have an instance of ITrackCancel but I can't work out how to show the default GP progress dialog from there.  I have this as code at the moment but it still doesn't show the dialog. Any ideas?

        IStepProgressor pStepPro = (IStepProgressor)TrackCancel;
        pStepPro.MinRange = 0;
        pStepPro.StepValue = 1;
        pStepPro.Message = "Starting Import...";
        pStepPro.Position = 0;
        pStepPro.Show();

        ...

        pStepPro.Hide();


texts: Font Formatting in Advanced Label Expression ArcGIS/Python

texts:
Almost have this quandry figured out, but inevitably a new problem has popped up. I finally have the "code" part down to do what I want, but now the font formatting isn't acting properly. When I use a simple expression with font tags I can create this image.

However, when I use the advanced expression as follows, it returns the garbled image after. I'm using the tags the same way and can't find any reason they shouldn't be working. Without quotes I get a syntax error. EDIT - Trimmed code to focus on issue area. Font tags were improperly placed, see answer 1.

if float([Percentage_Complete]) >= 90:
lab1 =  "<bol><clr green='255'>" + lab1 + "</bol></clr>\n"
if float([Percentage_Complete1]) >= 90:
lab2 =  "<bol><clr green='255'>" +  lab2 +  "</bol></clr>\n"
if float([Percentage_Complete2]) >= 90:
lab3 =  "<bol><clr green='255'>"+  lab3 + "</bol></clr>\n"
if float([Percentage_Complete3]) >= 90:
lab4 =  "<bol><clr green='255'>"+  lab4+  "</bol></clr>\n"


texts: Calculating area per range of values (classes) in raster using ArcGIS?

texts:
My data consists of several points each with xy coordinates and height.

I used spatial analyst tools -> interpolation -> natural neighbor to get an interpolation map for the height distribution, separated in categories like 10m-20m, 20m-30m and so on. This works fine.

Now, I would like to get the area of those categories. For example, I want to see that the interpolated area in the 10m-20m range is 200m².


texts: Writing If-Then Statement in ArcGIS Field Calculator using Python Parser?

texts:
Scenario:

I have four fields: A, B, C and D

Field A is text-based and contains values from '000' up to some number 'xxx'
Field B contains a set of values (also text, but represented by numbers)
Field C contains a set of values (also text, but represented by numbers)
Field D is empty and waiting to be populated

I would like to construct an if-then statement in Field Calculator that does something like this:

If Field A = 000 then Field D = Field B,
If Field A =/= 000 then Field D = Field C

I have spent the better part of my day trying to find the proper 'language' to ask Field Calculator, but I'm stuck.  I have 0 programming experience also.

Could anybody show me how this should be set up either in Python?

I'm using ArcGIS 10.0.


texts:
Is anyone working on porting the Canvec Symbology toolkit to ArcGIS 10? How would you go about it if you were?

The canvec symbology toolkit is comprised of a canvec.style file (which is an ms-access mdb with a different extension, containing tables of the style names and embedded symbology bitmaps), and an arcmap composition with some VBA code to loop through the table of contents and apply the styles as appropriate.


texts:
I am working with Digital Elevation Models in format TIFF 2km * 2km. As far as I know DEMs should have a 3d function but in my case the DEMs are flat. They only look 3d because of the hill shading. But when I open the DEM in ArcScene it is flat.

How can I turn the DEM into a real 3d model? How can I see the z values?


texts:
I wish to download the raster image of precipitation data.
Is that possible? I think I'm doing something wrong when I export it.
Here is my code:
//Rain
var lluvia = ee.ImageCollection('JAXA/GPM_L3/GSMaP/v6/operational')
.map(function(im){return im.clip(Uy);})
.filter(ee.Filter.date('2021-01-25', '2021-02-02'));
print(lluvia.size());
var precipitation = lluvia.select('hourlyPrecipRate');
var cumPrec = precipitation.sum();
Map.addLayer(cumPrec, vislluvia, 'Lluvias');

Export.image.toDrive({image: cumPrec,
description: "Imagen",
folder: "gee_data",
fileNamePrefix: "Imagen",
region: Uy,
scale: 30,
shardSize: 100,
fileDimensions: 5000,
fileFormat: "GeoTIFF"

})


texts: Rename postGIS table and updating the GeoServer?

texts:
I want to rename one of my postGIS tables. However, when I do this by applying a simple

ALTER TABLE old_name RENAME TO new_name;

The spatial data, as I understand it, isn't updated with the new name. First, can these references be fixed by running the function:

SELECT probe_geometry_columns()

Is the table itself OK after this? Or should other measures be taken also?
Second, how do I update GeoServer so that it finds this new table? As of now, if I run the above commands, and edit the layer info on the GeoServer to point to the new name... I still get an error (GeoServer still tries to find it under the old name). Exactly, where should the layer name be changed on the GeoServer?


texts: Using arcpy.SpatialJoin_analysis to only return 'Join_Count'

texts:
I want to join 2 shapefiles, and receive a new output feature that only contains the count in the joins.

I know to use arcpy.FieldMappings() to set what fields map from the target and join features, and by default a Join_Count field is added to the output feature.

If I supplied an empty arcpy.FieldMappings() object to the SpatialJoin_analysis() function parameter, would I be able to just do a join and only get the Join_Count field?


texts: Labeling Excluded Data of Classification in ArcGIS legend?

texts:
I'm symbolizing a feature class by graduated colors classified with 8 classes, but need to symbolize Excluded Data (Nulls for the classified field). However, the excluded data doesn't label, or at least follow the conventional labels I've set up based on SQL queries for a separate field.

Is there any way to label these excluded features based on a different field (not the null classified field) without having to duplicate the layer?


texts: QGIS crashes when trying to return layer from processing toolbox console

texts:
The following script is intended to take a polygon layer and extract the features having a particular value ('natural') in a given field ('landuse') and return the selection.

Everything seems to work fine until the final line where the results are returned. QGIS crashes.

If i set Output to "none", the script runs fine. And I check the resulting layer with vlayer.isValid and it always returns true.

Moreover, if instead of 'OUTPUT': 'memory:', i set 'OUTPUT': 'results.shp',  the resulting layer, written on disk, is good.

I don't know what's causing all the trouble. To recreate it you can pick any polygon layer (be it geopackage or shapefile) that has a selection attribute and set the corresponding values under 'FIELD' and 'VALUE'.

from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                   QgsProcessingAlgorithm,
                   QgsProcessingException,
                   QgsProcessingOutputNumber,
                   QgsProcessingParameterDistance,
                   QgsProcessingParameterFeatureSource,
                 QgsProcessingParameterVectorDestination,
                QgsProcessingParameterRasterDestination,
                QgsProcessingParameterFeatureSink ,
                QgsProcessingOutputVectorLayer,
                QgsVectorLayer
                   )
import processing

class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):
def tr(self, string):
    return QCoreApplication.translate('Processing', string)

def createInstance(self):
    return ExampleProcessingAlgorithm()

def name(self):
    return 'bufferrasterextend'

def displayName(self):
    return self.tr('Buffer and export to raster (extend)')

def group(self):
    return self.tr('Example scripts')

def groupId(self):
    return 'examplescripts'

def shortHelpString(self):
    return self.tr('Example algorithm short description')

def initAlgorithm(self, config=None):
    self.addParameter(
        QgsProcessingParameterFeatureSink(
            'OUTPUT',
            self.tr("Output layer"),
            QgsProcessing.TypeVectorAnyGeometry
        )
    )

def processAlgorithm(self, parameters, context, feedback):
    extract_natural = processing.run("native:extractbyattribute", {
        'FIELD': 'landuse',
        'INPUT': 'path_to_layer/vienna_landuse.gpkg|layername=vienna_landuse',
        'OPERATOR': 0,
        'OUTPUT': 'memory:',
        'VALUE': 'natural'
    })

    vlayer = extract_natural['OUTPUT']

    if vlayer.isValid():
        feedback.setProgressText("Layer ok")
    else:
        feedback.setProgressText("Layer is wrong!")

    return {
        #'OUTPUT': None # works fine
        'OUTPUT': vlayer # qgis crashes
    }


texts: Simple way to have an offline web based map?

texts:
My goal is to have a world map embedded in a web page, but without access to the internet. Preferrably only a set of tiles in combination with some slippy map js-lib. Storage space is not really an issue, and max zoom level is not required.

What is the easiest way to accomplish this?

I found a project that seems to fit the bill - Maps on a stick, but it seems to be dead?


texts: i have ab Excel table which it's fields are (year,month,day,hour,AM/Pm), how to convert those Fields into a Date format table in Arcgis and then creating a time frame Visualization.


texts:
I need to process raster data stored on HDFS in GeoTIFF format using GeoTrellis. Each file is composed by 8 bands and I need to perform NDVI and other map algebra operations selecting specific bands.

I'm able to handle input files with a MultibandGeoTiff object:

val geoTiff = HadoopGeoTiffReader.readMultiband(new Path("path/to/tiff/file"))(sc)

and then do some processing with the function tile.combineDouble(). Once I have a processed tile, how can I save it on an output file (in GeoTIFF format) on HDFS or on a specific host's filesystem?


texts: Exporting Layer to Feature Class is Incomplete without Error

texts:
I have a CSV of 1.2 million records that I want export to a GDB as a feature class. I've done this numerous times, but this was before upgrading from 10.5. Steps:

Add CSV
Display X/Y table
Export to GDB

The result (without error) is a FC in the GDB, but with 234,000 records instead of the 1.2 million in the CSV. I also tried exporting it off as a shapefile, knowing that there's a 2GB limit on those. Is there something I need to change in the ArcMap config files to increase the memory? I see several posts out there about changing the machine's virtual memory, but I've never had to do that for this.

I'm running ArcMap 10.6.1 on a 64-bit Windows 10 Pro machine.

EDIT: I have to import the CSV to a GDB and then display the X/Y in order to make it work. I used to be able to just import the CSV and display XY off the bat.


texts: In Vertex Tool (QGIS 3.2.2) How can I select the next (or previous) vertex?

texts:
I want to correct a duplicate vertex by deleting one.

However, if I use click and drag, I select both and delete both.  To select only one, I could select the previous (or next) vertex and navigate to the next (or previous) one and then delete.

I can't find a keyboard shortcut or menu option to navigate-select to next vertex.


texts:
I need to know how to detect which city has certain coordinates.
Which is the best way to do this?


texts: Georeference vertical images on 3D wireframe

texts:
I would like to make a 3D wireframe model of a room in a building and then georeference images of the walls on the wireframe. Once I have the images rendering on the walls I want to insert points on the walls that are linked to an Access Database containing a range of attributes.
Can anyone assist with the workflow please.


texts:
I have a set of regions and would like to quantify the degree to which the boundary of these regions are situated on/near rivers. I have turned the regions from polygons to polylines in arcgis 10.1 and have the rivers as polylines.

Is there any way to quantify what percentage of the boundaries fall on rivers, and therefore statistically show the degree of overlap between these two vector features?


texts: Creating polygonal windows for spatial analysis in r

texts:
I am attempting to normalize a spatial point pattern using a population density raster in r. My problem is that the function I am trying to run is recognizing the "white space" around the raster I am using as "NA", giving me the error message: Values of the covariate ‘Y’ were NA or undefined at 32% (19493 out of 60241) of the quadrature points. Occurred while executing: ppm.ppp(X, ~offset(log(Y)), covarariates = list(Y)).

Is there a way to create a window or data frame that traces the outline of my raster set? It should be noted that I am a grad student with very limited experience using r.

    library(spatstat)
    library(maptools)
    library(raster)
    library(geostatsp)
    library(rgdal)

    X <- readShapePoints("~/Documents/Data/Cyber/US/CONTIG_US/CONUS02_FREQ")
    Y <- raster("~/Documents/Data/Cyber/US/US_POP1.tif")

    ##Set the extent

    e <- extent(-125,-67,25, 50)

    ## create a window for spatstat

    WIN <- as(e,"SpatialPoints")
    X.ppp <- as(X, "ppp")

    ## convert raster to an im object for spatstat

    Y <- asImRaster(Y)
    Y <- eval.im((Y))

    Y <- as.im(Y)

    ## Unmark X for model:

    X<- unmark(X.ppp)

    Z <- ppm.ppp(X, ~offset(log(Y)), covarariates=list(Y))

The raster plot looks like this:


texts:
I'm trying to parallelize with multiprocessing.Pool() some processes that involve queries of R-trees. In the non-parallel procedure, it works as expected.

But when I try to do it in parallel (even with an example of running only 1 process in 1 core) the same query returns empty results.

I'm really puzzled about this, it goes beyond any logic for my understanding.

EDIT

Here is the simplest example that I could reproduce:

from multiprocess.pool import Pool
from rtree import index
class Test(object):
    def __init__(self):
        self.idx = index.Index()
        left, bottom, right, top = (0.0, 0.0, 10.0, 10.0)
        self.idx.insert(0, (left, bottom, right, top))

    def test(self, a):
        print(list(self.idx.intersection((1.0, 1.0, 2.0, 2.0))))

b = Test()
p = Pool()
res = p.map(b.test, range(3))
b.test(1)

output:

[]
[]
[]
[0]


texts:
I am new to arcpy. I am attempting to batch clip over 350 rasters to one the extent of one mask raster. I tried it in model builder and it wouldn't run so now I am  trying it in arcpy but I continue to get error messages.

Here is my code:

import arcpy
import os
import glob
from arcpy import env
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *

#Sets environment
env.workspace = "E:/AlexandraOrrego/GIS_2017/RasterData/Reclass2.gdb"

mask = "E:/AlexandraOrrego/GIS_2017/Shapefiles/westus_forest"
outputDir = "E:/AlexandraOrrego/GIS_2017/RasterData/clipped.gdb"

rasList = arcpy.ListRasters()
for raster in rasList:
    outRasterName = raster.replace("r2", "c")
    rasterObject = arcpy.sa.ExtractByMask(raster, mask, outRasterName)
#Saves Raster
    rasterObject.save(os.path.join(outputDir, outRasterName))
#prints once each raster layer has been clipped
    print outRasterName + " Clipped Successfully!"

 Runtime error
Traceback (most recent call last):
  File "<string>", line 17, in <module>
TypeError: ExtractByMask() takes exactly 2 arguments (3 given)
>>>


texts: How to undo a point while digitizing with QGIS 2.18? I am digitizing polygon features.


texts: FieldName underscore corrupts Label Expression for Font Size

texts:
I am trying to change the size of a secondary line in the Label Expression box:

  MID( [TRACT] ,4,11) & vbnewline & [FERC_NAME]

This approach does not work but it DOES work when the field name is something like TRACT.

  "" & [FERC_NAME] & ""

Any ideas as to what is going on?


texts: Merging contour lines generated by gridded data

texts:
I have a set of gridded data, which I have used to generate contour lines based on separate values in the data.

In the example provided the contour lines shown correlation coefficients.  The contours do a good job of interpolating between the grid points to indicate how good the correlation is to any of the three stations.

However, where contour lines overlap it is difficult to tell whether the correlation is stronger towards the blue or the red station. I would like to merge my contour lines together so that it is easy to tell which station has the best correlation for any given spot on the map.

Using the gridded data I sorted by the strongest correlation and created a convex hull around each as shown below.  However, this doesn't interpolate between the points like the contours do.

Is there a way to merge my contours together to create a smooth boundary and show where the highest correlation occurs for each?


texts:
I want to be able to select by attributes for dates that change weekly. I have gotten an invalid expression error for the follow script where the date needs to change to the current date and the same date just the following year every time I run it. I want to be able to create a model with multiple versions of this script.

import arcpy
in_layer = arcpy.GetParameterAsText(0)
in_field = arcpy.GetParameterAsText(1)
one_date = arcpy.GetParameterAsText(2)
two_date = arcpy.GetParameterAsText(3)

new_selection = "NEW_SELECTION"
where_clause = "in_field >= date 'one_date 00:00:00' AND in_field <= date 'two_date 00:00:00'"

arcpy.SelectLayerByAttribute_management(in_layer, new_selection, where_clause)


texts: Printing name of field with highest value?

texts:
I am trying to find the field with highest area value in each row of the 5 fields. Now I have found the highest value. But it doesn't really do me any good if I couldn't find which field is that max.

I have some rudimentary understanding and practice of Python. The very similar answer here:

need to find max value among fields AND PRINT FIELD NAME

When I applied the code, it always gave me a failure, which is probably due to my lack of understanding of Python.

  Error 000539, Python syntax.

The code I adapted from that question is shown as below:

import arcpy

table = "L:\......\ParcelSummary.dbf"
fields = [ !SUM_acres_! , !SUM_acres1! , !SUM_acre_1! , !SUM_acre_2! , !SUM_acre_3! ]
n = len(fields)

maxfield = "Highest"
maxname = "HighestF"
arcpy.AddField_management(table, maxfield, "FLOAT")
arcpy.AddField_management(table, maxname, "TEXT")

fields2 = fields[:] # shallow copy
fields2.extend([maxfield, maxname])

with arcpy.da.UpdateCursor(table, fields2) as cursor:
for row in cursor:

    check = row[:n]
    maxval = max(check)

    row[-2:] = maxval, fields[check.index(maxval)]

    cursor.updateRow(row)

Also, a pretty basic question about this: shall I put the code into the Field Calculator expression or shall I put it in Pre-Logic Script Code? If I put them in the Codeblock, then what shall my expression be?

As for the field names, the ones shown in the picture are the aliases of the fields instead of names.

My exact process and the error message I got.

Seems like the error is from line 1

The picture above shows the first few lines of the expression, what could have gone wrong here?


texts: ArcSDE : Cascading Subtypes in ArcGIS 10/SQL Server

texts:
I'm currently creating an ArcSDE DB on an SQL Server. It is designed to complete species survey, with information from the field.

In order to help themselves, they asked me to implement a particular yet logical behavior: in a first field Order, they want a list of related taxa. In a second field, a list of Genera (Genus) related to their Order choice. In a third field Species the list related to their Genus choice.

Basically, it is the subtype/domain behavior, but here, the field Genus should be a subtype and a domain at the same time and it's not possible.

So, I thought about implementing this behavior in SQL Server but is it possible? Moreover, will this behavior work directly in ArcMap when working in a editing session on the DB?

Other ideas?
Thanks, Merlin


texts: No layer available in the response of GetCapabilites request

texts:
I am doing a GetCapabilites request on geoserver to get the bounding box of wms layer but in the response, i am only getting a bunch of EPSG projections but no layers?
this is the request i am using:

http://localhost:8080/geoserver/wms?service=WMS&request=GetCapabilities&version=1.0.0

Also i have read the documentation of get capabilities https://docs.geoserver.org/latest/en/user/services/wms/reference.html#wms-getcap
and they have mentioned a namespace parameter to bound the result of GetCapabilites request to a specific layer,what should be the value of this parameter ?


texts:
I am trying to use data from GapMinder. I have downloaded the Excel spreadsheet data for CO2 emissions and for global GDP. I have uploaded the file to my folder connection in ArcGIS Pro, and dragged the table to the Contents pane and then joined the table to a World Countries layer from the AllPortal by "Country". When the join is complete, many countries are missing, like: Russia, Greenland, Venezuela, and the Democratic Republic of the Congo.

Has anyone ever tried using GapMinder data in ArcGIS Pro? If so, how were you successful in getting all of the countries to show on the map?


texts: Access attributes of a deleted IFeature in ArcObjects 10.1

texts:
I am migrating code for an ArcGIS extension from 9.3 to 10.1.
The code is throwing an exception in one place.

My duty is just to migrate it to 10.1, so I really don't like to do any modification more than few lines of code.

The code is searching for a feature that was deleted by the user on a child version,
but
1 - it is not posted yet to the default version
2 - The edit session is still open and not saved.

so before the save it is trying to validate some business rules.

IFeatureWorkspace childWorkspace = (IFeatureWorkspace)childVersion;
IFeatureClass childOperationalUnit = childWorkspace.OpenFeatureClass("OperationUnit");
IFeature childFeature = childOperationalUnit.GetFeature(Id_Of_Deleted_Feature);
if (childFeature != null)
{
 string operationName = childFeature.get_value(Index_Of_Name_Column);

if (operationName == "somevalue")
......
}

So, the code is trying to access the attributes of a deleted feature.
The code was working on 9.3

But on 10.1 is not working and giving me the error:

System.Runtime.InteropServices.COMException (0x8004100C): The feature is deleted
The error happening on the line where we try to get the attributes of the feature with get_value

so, 10.1 can find the feature but cannot access its attributes.

is there any easy way to keep the code working without really modifying lots of code?


texts:
I downloaded a tool and installed it in the ArcGIS, which can allow me split by attribute very easily and quickly. I wonder how I can use it as the built-in tools of the ArcGIS, so that I can write python code to use it as I need? For example, we can use Clip as writing arcpy.Clip_Managerment. I want to have a similar way of writing code for the split by attribute tool installed in the ArcGIS.

I am using ArcGIS 10.3.1 Advanced License.
The tool website is http://www.umesc.usgs.gov/management/dss/split_by_attribute_tool.html.
The following picture shows the tool I installed in the ArcGIS.


texts: Sql to view points of intersection with line not working

texts:
I have two tables one called test lines and other points called testepoint with some geometries drawn in QGIS to test my sql. Are as in the image below:

I want to select the points that inteceptam a row, the result should return me the id of these three points:

This is my sql that returns no records:

    Select testepoint."id"
    FROM
    testepoint
inner JOIN
    teste
    ON ST_Intersects(teste.geom, testepoint.geom)

Can anyone help me on what I might be doing wrong?
I drew the shapes in QGIS using Google Streets layer based on EPSG: 3857 and carried in PostGIS

I ran the sql:

    Select testepoint."id"
FROM
testepoint
inner JOIN
teste
ON ST_DWithin(teste.geom, testepoint.geom, 1)
WHERE
teste."id"= 1

The result:

But I need this order: 1,2,4
The order of alignment.


texts: If I take a Natural Earth shapefile and delete all countries except one - This will affect that original shapefile. Can I take that shapefile duplicate it and use that shapefile to delete all countries except one and then save that without affecting the original shapefile?


texts: Displaying of bands for multiple datasets in Google Earth Engine

texts:
I am classifying landsat 8 and sentinel 1 in google earth engine but the console tells me i do not have the required bands for classification even after assigning all the bands.
Please find attached a link to my work flow
......... javascript for GEE
var season1 = Landsat.filterDate('2018-10-01', '2018-12-31')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season1:', season1);

var season2 = Landsat.filterDate('2019-01-01', '2019-03-31')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season2:', season2);
var season3 = Landsat.filterDate('2019-04-01', '2019-06-30')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season3:', season3);
var season4 = Landsat.filterDate('2019-07-01', '2019-09-30')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season4:', season4);

// Define False color visualization
var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B5', 'B4', 'B3'],
  gamma: 1.7
};

///Visualize the FCC of each season Landsat images
Map.addLayer(season1, rgbVis, 'Landsat MSI - season1',false);
Map.addLayer(season2, rgbVis, 'Landsat MSI - season2',false);
Map.addLayer(season3, rgbVis, 'Landsat MSI - season3',false);
Map.addLayer(season4, rgbVis, 'Landsat MSI - season4',true);

// Below four lines select season wise Sentinel 1 RADAR data again median and clipped to MIS
var vh_s1 =  S1.filterBounds(mis).filterDate('2018-10-01', '2018-12-31').median().clip(mis).select('VH').rename(['vh_s1']);
var vh_s2 =  S1.filterBounds(mis).filterDate('2019-01-01', '2019-03-31').median().clip(mis).select('VH').rename(['vh_s2']);
var vh_s3 =  S1.filterBounds(mis).filterDate('2019-04-01', '2019-06-30').median().clip(mis).select('VH').rename(['vh_s3']);
var vh_s4 =  S1.filterBounds(mis).filterDate('2019-07-01', '2019-09-30').median().clip(mis).select('VH').rename(['vh_s4']);

// Here we combine all the seasonal Landsat and Sentinel 1 to a single image collection called
// 'dataset'
var dataset = season1.addBands(season2).addBands(season3).addBands(season4).addBands(vh_s1).addBands(vh_s2).addBands(vh_s3).addBands(vh_s4);
print('dataset:', dataset);

//Below codes split the training data into Training (70%) and Validation (30%) - Random
var trdata_table = field.randomColumn({seed: 1});
var training = trdata_table.filter(ee.Filter.lt('random', 0.70));
var validation = trdata_table.filter(ee.Filter.gte('random', 0.70));

// Create a new feature collection (vector) with all the bands in 'dataset' sampled to training points
var training = dataset.sampleRegions({
  collection: training,
  properties: ['class_id'],
  scale: 20,
  tileScale: 8
});
print('training:', training);
// Create a new feature collection (vector) with all the bands in 'dataset' sampled to validation points
var validation = dataset.sampleRegions({
  collection: validation,
  properties: ['class_id'],
  scale: 20,
  tileScale: 8
});

//Define bands require for training and classification. So Bands 1,2,3,4,5,6,7,9,10,11 for 4 seasonal
// images and Sentinel 1 (1 band) representing 4 seasons
var bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B9', 'B10', 'B11', 'B1_1', 'B2_1', 'B3_1', 'B4_1', 'B5_1', 'B6_1', 'B7_1', 'B8_1', 'B9_1', 'B10_1', 'B11_1', 'B1_2', 'B2_2', 'B3_2', 'B4_2', 'B5_2', 'B6_2', 'B7_2', 'B8_2', 'B9_2', 'B10_2', 'B11_2', 'B1_3', 'B2_3', 'B3_3', 'B4_3', 'B5_3', 'B6_3', 'B7_3', 'B8_3', 'B9_3', 'B10_3', 'B11_3', 'vh_s1', 'vh_s2', 'vh_s3', 'vh_s4'];
print('bands:', bands);


texts: Mapbox Vector Tiles in Leaflet - How to get feature info?

texts:
I'm using geoserver to serve PostGIS layer (contains points only) as MVT (Mapbox Vector Tiles). I have a simple Leaflet application that displays that layer along with OSM or Google satellite imagery background.

Here is the code:

var osmUrl = 'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
osmAttrib = '&copy; OpenStreetMap contributors',
osm = L.tileLayer(osmUrl, { maxZoom: 18, attribution: osmAttrib }),
map = new L.Map('map', { center: new L.LatLng(42.647, 18.101), zoom: 13 }),
drawnItems = L.featureGroup().addTo(map);

// MVT layer style
var vectorTileStyling = {
    trgovinas: {
                icon: new L.Icon({
                                    iconUrl: 'icons/shop.png',
                                    iconSize: [20, 20], // size of the icon
                                    iconAnchor: [10, 10], // point of the icon which will correspond to marker's location
                                })
                }
    };

// MVT layer url
var trgovineUrl = "localhost:8080/geoserver/gwc/service/tms/1.0.0/volonteri:trgovinas@EPSG%3A3857@pbf/{z}/{x}/{-y}.pbf";

// MVT layer options
var mapboxVectorTileOptions = {
    rendererFactory: L.svg.tile,
    attribution: '&copy; OpenStreetMap contributors, &copy; MapBox',
    vectorTileLayerStyles: vectorTileStyling,
};

var trgovine = L.vectorGrid.protobuf(trgovineUrl, mapboxVectorTileOptions);

// Layer control + layers definition
L.control.layers({
    'osm': osm.addTo(map),
    "google": L.tileLayer('http://www.google.cn/maps/vt?lyrs=s@189&gl=cn&x={x}&y={y}&z={z}', {
        attribution: 'google'
    })
}, { 'drawlayer': drawnItems ,
     'trgovine': trgovine.addTo(map)}).addTo(map);

What I want to do is to store features (geometry and attributes) to a variable when I click a point on map so I can make a popup that contains information about that point (from MVT layer).

This should be simple because data is already on client side (that is why I'm using MVT), but I'm new to Leaflet and JavaScript events.


texts: How to write a custom map canvas item by subclassing QObject and QgsMapCanvasItem?

texts:
I'm writing a QGIS python plugin where I place some QgsVertexMarkers on the canvas like this:

marker = QgsVertexMarker(self.canvas)
marker.setCenter(QgsPoint(point[0],point[1]))

Now i would like being able to drag and drop them. I tried implementing a class that inherits QObject and QgsMapCanvasItem so I can use mouse events.

class myVertexMarker(QObject,QgsMapCanvasItem):
    def __init__(self, canvas):
        QObject.__init__(self)
        QgsMapCanvasItem.__init__(self,canvas)
        ...

    def setCenter(self,point):
        self.center = point
        pt = QPointF(self.toCanvasCoordinates(self.center))
        self.setPos(pt)

    def some setters, mouse events, paint and so on...

Now when I create an object of that class and call setCenter I get this error message:

pt = QPointF(self.toCanvasCoordinates(self.center)) #TypeError: could not convert
'myVertexMarker' to 'QgsMapCanvasItem'

What am I doing wrong?


texts:
I am opening an Excel spreadsheet in MapInfo.

One column of the Excel spreadsheet holds a date per row. Some of the rows dont have an assigned date this is fine.

When opening the table in MapInfo, the browser window shows that these blank cells have now been infilled with the date 1899?

How can I stop this happening?


texts: I obtained some elevation data from the USGS website and tried importing it into BlenderGIS and SVDesigner (SoilVision). However, I get a weird render

texts: The first image is rendered on SVDesigner, the second image is rendered on Blender GIS and the third image is the raster file I obtained from the USGS website.


texts: How to generate a thumbnail of an openlayers map?

texts:
I would like to allow the user to generate a thumbnail image of what the user is seeing in the OpenLayers map canvas so the user may print, save, share, email, or whatever.

Any ideas of how i would go about doing this?

I'm looking for solutions that generate an image. I don't want to save the innerHTML, feature data or anything other than an image. An image is much easier to share and cannot later be modified.

Examples:

or this will also work:


texts: I have a KML layer of bicycle routes and a street line layer.  Using either ArcGIS or QGIS, how can I select the streets that correspond to the bicycle routes?


texts: Get coordinate from mouse click in ArcMap using ArcObjects?

texts:
I am using ArcObjects C# tool in ArcGIS that identifies map coordinates when a user clicks on a location within a map. Here's the code that I have used:

My map use units of Meters so the tool displays the Meters coordinates. How can I have this tool return degrees decimal?


texts: QGIS loads PostGIS layers but DBManager can't connect

texts:
DB Manager can't connect to PostgreSQL, but "Add PostGIS Layer" can.

My eventual goal is to save created layers to PostGIS. I assume I can do this using the DBManager Tool. I can alternatively use ogr2ogr or shp2pgsql, but I was hoping QGIS would make workflows between shapefiles and PostGIS a bit easier.

This behavior seems odd. Why can it connect one way and not the other?

Side Note: Why does DBManager display my password back to me?

I can connect easily with psql on command line

Additional Info:

This is a local database
When I originally connected using 'Add PostGIS Layer', I chose not to save the username and password. I am connecting using my OS user account.
I'm running OS X 10.9.2, PostGIS 2.1, PostgreSQL 9.3.3, and QGIS 2.2.0
PostgreSQL, PostGIS and QGIS were all installed using Homebrew. See Homebrew install info below.

Homebrew info:

brew info qgis-22
qgis-22: stable 2.2.0, HEAD
http://www.qgis.org
/usr/local/Cellar/qgis-22/2.2.0 (14074 files, 493M) *
  Built from source with: --with-grass, --with-api-docs, --with-postgis
From: https://github.com/dakcarto/homebrew-osgeo4mac/commits/master/Formula/qgis-22.rb
==> Dependencies
Build: cmake ✔, bison ✔
Required: qt ✔, pyqt ✔, qscintilla2 ✔, qwt ✔, qwtpolar ✔, gsl ✔, sqlite ✔, expat ✔, proj ✔, spatialindex ✔, fcgi ✔, gdal ✔, pyspatialite ✔
Recommended: postgresql ✔, gpsbabel ✔
Optional: postgis ✔, grass ✔, qt-mysql ✘, orfeo ✘, r ✘, saga-gis ✘

brew info postgis
postgis: stable 2.1.1, HEAD
http://postgis.net
/usr/local/Cellar/postgis/2.1.1 (44 files, 8.3M) *
  Built from source
From: https://github.com/Homebrew/homebrew/commits/master/Library/Formula/postgis.rb
==> Dependencies
Build: autoconf ✔, automake ✔, libtool ✔, gpp ✔
Required: postgresql ✔, proj ✔, geos ✔, json-c ✔, gdal ✔

brew info postgresql
postgresql: stable 9.3.3 (bottled)
http://www.postgresql.org/
Conflicts with: postgres-xc
/usr/local/Cellar/postgresql/9.3.2 (2968 files, 66M)
  Poured from bottle
/usr/local/Cellar/postgresql/9.3.3 (2967 files, 65M) *
  Poured from bottle
From: https://github.com/Homebrew/homebrew/commits/master/Library/Formula/postgresql.rb
==> Dependencies
Required: readline ✔
Recommended: ossp-uuid ✔


texts: GeoServer 2.10.1 Center Justify Polygon Label Style File .xml

texts:
The Title kind of says it all.  Trying to Center Justify a label of a polygon using a sld style file on Geoserver 2.10.1.  Here is the code snippet followed by a screenshot of the label now.

          <sld:Rule>
        <ogc:Filter>
          <ogc:PropertyIsLike wildCard="*" singleChar="#" escapeChar="!">
            <ogc:PropertyName>building</ogc:PropertyName>
            <ogc:Literal>*</ogc:Literal>
          </ogc:PropertyIsLike>
        </ogc:Filter>
        <sld:MinScaleDenominator>0</sld:MinScaleDenominator>
        <sld:MaxScaleDenominator>3600</sld:MaxScaleDenominator>
        <sld:Name>Yes</sld:Name>
        <sld:Title>Building</sld:Title>
        <sld:PolygonSymbolizer>
          <sld:Fill>
            <sld:CssParameter name="fill">#d4d6e4</sld:CssParameter>
          </sld:Fill>
          <sld:Stroke>
            <sld:CssParameter name="stroke">#d4d6e4</sld:CssParameter>
            <sld:CssParameter name="stroke-width">0.5</sld:CssParameter>
            <sld:CssParameter name="stroke-linejoin">round</sld:CssParameter>
          </sld:Stroke>
        </sld:PolygonSymbolizer>
        <sld:TextSymbolizer>
          <sld:Label>
            <ogc:PropertyName>name</ogc:PropertyName>
          </sld:Label>
          <sld:Font>
            <sld:CssParameter name="font-size">12</sld:CssParameter>
            <sld:CssParameter name="font-family">Nimbus Sans L Regular Condensed</sld:CssParameter>
          </sld:Font>
          <sld:LabelPlacement>
            <sld:LinePlacement>
              <sld:PerpendicularOffset>0</sld:PerpendicularOffset>
            </sld:LinePlacement>
          </sld:LabelPlacement>
          <sld:Fill>
            <sld:CssParameter name="fill">#749dae</sld:CssParameter>
          </sld:Fill>
          <sld:VendorOption name="maxDisplacement">50</sld:VendorOption>
          <sld:VendorOption name="labelAllGroup">TRUE</sld:VendorOption>
          <sld:VendorOption name="removeOverlaps">TRUE</sld:VendorOption>
          <sld:VendorOption name="autoWrap">40</sld:VendorOption>
          <sld:VendorOption name="group">TRUE</sld:VendorOption>
          <sld:Halo>
            <sld:Radius>
              <ogc:Literal>0.26</ogc:Literal>
            </sld:Radius>
            <sld:Fill>
              <sld:CssParameter name="fill">#bce6f4</sld:CssParameter>
            </sld:Fill>
          </sld:Halo>
        </sld:TextSymbolizer>
      </sld:Rule>


texts: Is it possible to implement interactive labels in World Wind SDK. In API there is no labeling management interface for objects. Is there any usefull library that adds labeling mechanism to World Wind SDK.


texts: How to put value from different fields using Python?

texts:
Here I added a new field name "visArea" and want to put in the number from the "Count" column which "Value" equals one.That means the number of visible point(value=1)from selected observe spot.I used embedded Python 2.7 within ArcGIS 10.4 and my code are:

fc="PythonVisi2.tif"
Cursor=arcpy.da.UpdateCursor(fc,["Value","Count","visArea"])
for row in Cursor:
    if row[0] != 0:
    row[2]=row[1]
    cs3.updateRow(row[2])

But I got "TypeError: argument must be sequence of values".Which part goes wrong or how shold I fix the problem?


texts: Changing SRID of existing data in PostGIS?

texts:
While importing my shapefile data to PostGIS, I did not select the proper Projection.

How do I now change the SRID of the data, without transforming the Coordinates?


texts: Making species map using km grid in ArcGIS Desktop?

texts:
In ArcGIS Desktop 10.5, I am trying to make a species map to show the locations certain species are observed. I want to do this using kmgrid. I got quite far I think but now I am stuck.

I made a km grid for the netherlands giving every square kilometre a unique code. I also used this code in the database with the species. I couldn't join the two tables because in the species table te unique code was used multiple times so I made a pivot table to make sure all codes where unique.

After this I could link the data but now arc gis shows all species at once.

I can manualy select every species seperately but this will take a lot of time.

Is there another step I can take so I can use data driven pages for example?


texts:
I parse a geoJson file resulting in a list of osgeo.ogr.geometry-classes (all Polygons). I also have a new one which I want to dissolve with the rest.

I'm using the ogr package in python, it does something like Union, but I was unable to get satisfying results. Does anyone have a solution for dissolving the two polygons?

Edit (example)
A sample was asked in the comments, so here it is;
The sample-WKT is taken from Wicket. Then I simply do a union, with my expected results being one polygon representing the outline (since they are overlapping). Maybe I totally misunderstand the Union method, but would like to have a 'dissolved' result.

import ogr

# exported from; https://arthur-e.github.io/Wicket/sandbox-gmaps3.html
wkt_one = "POLYGON((109.86328125 -0.7745125438579877,107.40234375 -44.13885576756879,171.38671875 -47.327653995607086,165.05859375 10.071628685919348,109.86328125 -0.7745125438579877))"
wkt_two = "POLYGON((90.52734375 15.554252501699544,-157.32421875 -6.7355309886055625,-146.07421875 34.24813554589752,104.23828125 44.789632547614076,90.52734375 15.554252501699544))"

polyOne = ogr.CreateGeometryFromWkt(wkt_one)
polyTwo = ogr.CreateGeometryFromWkt(wkt_two)

print polyOne
print polyTwo
print "---"
print polyOne.Union(polyTwo)


texts: Geodjango transform introduces large errors

texts:
Apologies for the possibly naive question.  I'm trying to understand how Geodjango (running on PostGIS) deals with coordinate transformations, and am experiencing some unusual behaviour.

Let's say I have a very simple model with a geometry field, with SRID 27700 (OSGB):

from django.contrib.gis.db import models
class MyGeoClass(models.Model):
    point = models.PointField(srid=27700)

I create an entry in this coordinate ref system:

from django.contrib.gis import geos
m = MyGeoClass(point=geos.Point([527860, 184660], srid=27700))
m.save()

Now, if I load the entry and transform it to a different coordinate system and back again, it introduces a slight error

m = MyGeoClass.objects.all()[0]
m.point.coords
>>> (527860.0, 184660.0)
m.point.transform(4326)
m.point.transform(27700)
m.point.coords
>>> (527860.0007865971, 184659.99966225715)

No surprise there - I guess this is just numerical rounding error.  But if I instead leave the point in SRID 4326 and issue a save() call, Geodjango automatically transforms it for me, with a much larger error:

m = MyGeoClass.objects.all()[0]
m.point.coords
>>> (527860.0, 184660.0)
m.point.transform(4326)
m.save()
m = MyGeoClass.objects.all()[0]
m.point.srid
>>> 27700
m.point.coords
>>> (527861.48656512, 184654.1590452015)

So, the transform-transform pathway introduces an error of around 0.001 metres, whereas the transform-save pathway introduces an error of 1.4 metres in the x-dimension and 5 metres in the y dimension.


texts: Why use QgsFeatureRequest() when iterating?

texts:
I am unclear on why one should use QgsFeatureRequest() when iterating over a file’s features.

Let’s say I would like to print out all IDs in my currently active dataset.
Now, most examples will provide something like this:

lyr = iface.activeLayer()
prov = lyr.dataProvider()
for feature in prov.getFeatures(QgsFeatureRequest()):
    print feature.id()

Whereas the following, without referring to QgsFeatureRequest(), would yield the same result:

lyr = iface.activeLayer()
prov = lyr.dataProvider()
for feature in prov.getFeatures():
    print feature.id()

I’d assume that the right way is using QgsFeatureRequest(), as most quality sources include it, but can someone explain what it actually does, and why one should use it?


texts: Is there a way in QGIS 3 to get a selected feature WKT?

texts:
In QGIS 2.x, I used the getWkt  plugin: https://github.com/hugoledoux/QGIS-getWKT to get a WKT of a selected feature.

But that plugin isn't present in QGIS 3.x.  Was that capability built into core QGIS 3?  Or is this a case where we are waiting for the plugin to catch up?


texts: Different legend properties with ArcGIS API for JavaScript

texts:
I have a map template that I use to display a variety of feature layers and dynamic layers.
I have found a snippet of code to add the legend to my map:
//add the legend
dojo.connect(myMap, 'onLayersAddResult', function (results) {
    var layerInfo = dojo.map(results, function (layer, index) {
        return { layer: layer.layer, title: layer.layer.name };
    });

    if (layerInfo.length > 0) {
        var legendDijit = new esri.dijit.Legend({
            map: myMap,
            layerInfos: layerInfo},
            "legendDiv");
            legendDijit.startup();
        }
    });

The issue is, I can't figure out how to configure the names of the layers, labels, descriptions, etc. Most of the layers in my service are symbolized by Category, so, I want to be able to display those category names in the Legend.
This is what I want to show..or something like it, but with a category or layer heading

This is what shows up. No heading, just the categorized layer without a title or anything.

Any help would be appreciated.
Thanks


texts:
Suppose I have geographic coordinates of "Saratoga, California, USA" as

Latitude:   37°15.8298′ N
Longitude: 122° 1.3806′ W

I know from here that in case of latitude 1° ≈ 69 miles and that longitude varies:

1° longitude = cosine (latitude) * length of degree (miles) at equator.

How many miles is 1° longitude at longitude: 122°1.3806′ W?


texts: How to calculate the 3D length of line segment in QGIS?

texts: I've been looking through all of the GRASS commands in the QGIS 2.10 processing toolbox and it seems like one of these must be able to do this but I haven't found one yet. There is 'v.sample' which can sample a DEM at a point location. Is there a way to do something similar to sample a DEM when calculating the length of a line segment?


texts: Seeking open source route planning software?

texts:
Let's assume for this post that I know the amount of time I'll be spending at each location; the only unknown is travel time between locations.

I've been loading the addresses into Google maps, eyeballing the best route and rearranging the multiple endpoints, then copying travel times back into an Excel spreadsheet that calculates a timed itinerary for my day. Multiple iterations of a route can take 30 to 45 minutes to plan out.

From my research, I've found route planning and fleet management systems designed for delivery companies that are proprietary and cost-prohibitive. I've also found add-ins for GPSs that can plan routing for "up to five!" destinations. Assuming I have that particular GPS brand and pay for annual software updates.

I'm looking for something that would take a list of addresses, load it into Google maps, and spit back the most efficient route with addresses and travel times preferably as an Excel or CSV file. A plus would be the ability to specify "most efficient" by time or mileage, or a fuzzy-logic magic blend of the two. Open source is preferred.

Has anyone ever run into anything that might meet my needs?


texts: stall-on-demand for contraction hierarchy routing: explanation with sample code needed

texts:
I've implemented a contraction-hierarchy routing system, which now works well, and I now need to add the stall-on-demand technique to speed it up; as explained here it can reduce the search space dramatically: over 25-fold in the example given.

However, I'm finding it hard to understand the few explanations I have found on the web, and the sample code given in the cited document lacks context. I'd be grateful if anyone can point me to fuller and more carefully explained examples and code.

So far I've found some code which is reasonably easy to understand in Christian Vetter's MoNav project. So that partially answers my own question. But more information, and a step by step explanation, would still be useful.


texts: Configuring columns when adding CSV file to QGIS

texts:
I want to load a csv file into QGIS, and the data in QGIS as point layer. I have tried to use this function:

Layer ► Add layer ► Add Delimited Text Layer.

When Im loading the file name, the data comes like this in the box:

All the data is just in one column. Latude and longitude should be x,y values. I had GHI, Latitude and Longitude in different columns in excel, like this:

To make the csv file I saved as "CSV (comma delimited).

So how can I load the csv file correctly? I havent done a operation like this before.


texts: Extract linear features (roads) from DEM

texts: I have a 1x1meter DEM and the original LiDAR point cloud. I need to create a vector feature with roads. Are there any tools in ArcGIS to do this? Or any ideas for the algorithms to try to do this in R or Matlab?


texts: Using buffer from IDW in ArcGIS 3D Analyst for flood spread?

texts:
I have a line in the centre of a road, with the roads elevation attached to its properties.
By using Generate Point along lines --> Add Surface Information --> Adding 10 cm to that, and then IDW 3d analyst I've made a new spread. Shown on the picture below. Then I've done a CutFill tool (or alternatively you can you use Raster Calculator) to get where the areas will be flooded.

But is is possible to get it to take into account the existing elevation model, so it stops the spread/buffer and doesn't create water behind a barrier?

See below images for clarification. The ones with red circles around them are examples of areas that shouldn't be flooded. Is there a way these can be removed/filtered, my own idea is to somehow check if they have a connection to the centreline of the road.

I have access to Spatial Analyst, 3D Analyst, Arc Hydro and ArcMap with an Advanced level license. I have to do this for a large area, not just this test area that is shown here.


texts: Restricting photo size for attributes displayed in qgis2web?

texts:
I have created field for objects in a layer that allow me to input links that will display a photo on mouse-over for that object in a webmap made with the qgis2web plugin for QGIS v3.2.

How do I restrict the image size to some preset boundary? Currently, each photo is displayed at full extent. I don't want to resize the image to explicit dimensions, because some are oriented vertically, while others are oriented horizontally. I would like to preserve the original aspect ratio of the photos, while reducing their overall size.

Added 9/24/18: I am exporting as OpenLayers.


texts: Change marker icon on click using leaflet

texts:
I have a map with many (>100) markers on it.  I would like to make it so that clicking on these changes the icon to a highlighted version.  I've created two custom icons, one regular, and one highlighted.  I've gotten this to work with single markers, but can't find any way to set it up so each icon can be changed by clicking on it.

Here is working code for one icon:

var testmarker = L.marker([44.0, -73.0], {icon: unselectedIcon});

testmarker.on('click', function(e) {
  testmarker.setIcon(selectedIcon);
});

However I have many markers and use a for loop to place them:

//points is a 2 dim array with lat/long pairs
for (i = 0; i <= points.length-1; i++) {
    var marker = L.marker([points[i][0], points[i][1]], {icon: unselectedIcon});
}

marker.on('click', function(e) {
  marker.setIcon(selectedIcon);
});

I tried giving each marker a unique variable name by using an array:

for (i = 0; i <= points.length-1; i++) {
    testmarker[i] = L.marker([points[i][0], points[i][1]], {icon: unselectedIcon});
    testmarker[i].on('click', function(e) {
alert(e.latlng); e.setIcon(selectedIcon);
});
}

Which just gives me "undefined" in the alert.

Changing testmarker[i].on to this.map.on only gives an alert on clicking the map (no icon change).

Changing that entire line to:

testmarker[i].on('click', function(e) {
alert(e.latlng); testmarker[i].setIcon(selectedIcon);
});

Changes the last one (when any is clicked), while also giving an undefined alert.

How can I set the markers up so that each can individually be changed with a click, but without using hundreds of repeated code blocks for each one?


texts: OpenLayers/GeoServer layer zooms in and then out again, why?

texts:
I have a strange ZOOM behavior.  I am using GeoServer 2.6.1/OpenLayers 2.13.1.  I have a layer I created from a MrSid file and am serving as a wms layer.

Lets call it the "problematic" layer.  When I serve this "problematic" layer alone, everything is ok.  When I serve it with a base layer, I get a strange effect:  when I zoom in, the "problematic" layer zooms in, then zooms back out again.

If I click on a check box to hide the "problematic" layer (and only show the base layer which is at correct zoom) and then I click on the check box to show the "problematic" layer again, it aligns correctly with the base layer.

On the GeoServer side, I am "forcing declared" on both the base layer and the "problematic" layer.  The Native SRS and the Declared SRS are both set to EPSG:2236. On the OpenLayers side, I am using the displayProjection setting and have it also set to 2236 as recommend in other threads.

What could be the issue?

OpenLayers Code:

var map, baseLayer, problematicLayer;

        // pink tile avoidance
        OpenLayers.IMAGE_RELOAD_ATTEMPTS = 5;
        // make OL compute scale according to WMS spec
        OpenLayers.DOTS_PER_INCH = 25.4 / 0.28;

        function init() {

            // Setup Projection
            epsg2236 = new OpenLayers.Projection("EPSG:2236");

            // Setup Format
            format = 'image/png';

            // Setup Extent
            var bounds = new OpenLayers.Bounds(
                    858665.911628, 516256.950565,
                    870969.136, 527480.405429
                );

            // Setup Map Options
            var options = {
                maxExtent: bounds,
                maxResolution: 128.890625,
                projection: epsg2236,
                displayProjection: epsg2236,
                units: 'm'
            };

            // Setup Map
            map = new OpenLayers.Map('map', options);

            // Setup tiled base layer
            baseLayer = new OpenLayers.Layer.WMS(
                "sf:base_lines - Tiled",
                "http://server01:8085/geoserver/sf/wms",
                {
                    "LAYERS": 'sf:base_lines',
                    "STYLES": '',
                    TRANSPARENT: 'true',
                    format: format
                },
                {
                    visibility: true,
                    buffer: 0,
                    displayOutsideMaxExtent: true,
                    isBaseLayer: true,
                    yx: { 'EPSG:2236': false }
                }
            );

            // The Problematic Layer
            problematicLayer = new OpenLayers.Layer.WMS(
                    "sf:SID_2014 - Tiled",
                    "http://server01:8085/geoserver/sf/wms",
                {
                    "LAYERS": 'sf:SID_2014',
                    "STYLES": '',
                    TRANSPARENT: 'true',
                    format: format
                },
                {
                    visibility: false,
                    buffer: 0,
                    displayOutsideMaxExtent: true,
                    isBaseLayer: true,
                    yx: { 'EPSG:2236': false }
                }
            );

 map.addLayers([baseLayer, problematicLayer]);

// Add MousePosition Control
            map.addControl(new OpenLayers.Control.MousePosition());

            // Zoom to Max Extent
            if (!map.getCenter()) {
                map.zoomToMaxExtent();
            }
 }

Update #1

I just uploaded a video of the issue so you can get an exact idea of what I mean by "zooms in and then out again":

https://vimeo.com/119384955


texts:
I want to create a number of featureclass using a template featureclass

import arcpy

from arcpy import env
env.workspace = "C:\new_docs\BIO-FILES\GIS_prog\Lesson1"
arcpy.CreateFeatures_management('C:\new_docs\BIO-FILES\GIS_prog\results","Precip2009Readings.shp","POINT","Precip2008Readings.shp","Disabled","Disabled")

When I run the above code I get:

  Parsing error SyntaxError: EOL while scanning string literal (line 5)


texts: Adding KML linestring intersection points to relevant linestring coordinates using QGIS

texts:
I have a KML document with many LineStrings which I'm parsing and building a graph from in an application. Unfortunately not all LineStrings line up correctly (having the same vertex/coordinate), and I have to find the intersections and insert them in the correct place in both LineStrings' coordinates-array, as in the picture:

I've found all the intersection points by using QGIS -> Vector -> Analysis Tools -> Line intersections, but I'm stuck on the last step of inserting the intersection points into the two relevant LineStrings' coordinate-arrays, and I was hoping to be able to do this still in QGIS so that I can export it as a new KML document with the new points correctly included.


texts: How to improve the color matching between a scanned image of a map and the original?

texts:
When I scan a section of paper map using a Canon Lide 200 A4 scanner, the resulting image has colours which are not true to the original.
The most annoying is where the predominantly white paper background of a map is rendered as a greyish blue colour instead of pure white. Processing in ESRI ArcEditor 9.2 such as clipping the image so that it can be butt joined and geo-referenced, results in even darker shades, even less true to the original.

I need a very simple solution to this problem as I'm too old to start learning how to use Photoshop or similar image processing software. In finding a solution cost isn't the issue - The learning curve is.
Can anyone offer advice?


texts: Understanding flow directions with Indeterminate values on Geometric Network using ArcMap

texts:
I am planning a water distribution network on ArcMap. The water is coming from source and going to customers like following.
The green squares are customer service points. I write labels of AncllayRole fields on network. When I click the Set Flow Direction button on toolbox, arrows are appearing until we come to the junctions. But then the arrows does not appears so Indeterminate points are signed on pipe lines.
Why is that, and how can I solve this issue?


texts: Grass GIS table creation error during r.to.vect

texts:
When I execute r.to.vect the following error appears:

r.to.vect input=sreviLC815109_-120_-120 output=temp feature=area --o --q
DBMI-DBF driver error:

ERROR: Unable to create table: create table temp ( cat integer, value
   double precision)

Any ideas what this could be?
My main issue is that this error appears absolutely sporadically in a script while the action is performed on many files. I could shift the error message a bit to a later point in time by deleting the temp vector map before, but it still appears. I am working with Grass 6.4.3.

EDIT:
When I do the action manually, everything works fine.


texts: Creating outer glow effect depending on the graduated symbology color of a shapefile - QGIS

texts:
I am working with a line shapefile that represents the road network of a city. I used a graduated symbology to draw each roads depending on their attribute value.

Now, I want to get the outer glow of each lines for visualisation purposes.
Going to Symbology > symbol > effect only allows me to add a outer glow effect to the whole network at one with a definite color. I, however, would like to have the glow color depending on the line it is linked to.

Here's an example in case I was not clear enough:

This is the network colored depending on the roads attributes values:

And this is what I get when I add the outer glow effect to it which is not what I am looking for:

What I want is for the red roads glow to be red, for the blue ones to be blue, for the green ones to be also green etc.

Is there a way of doing so? I have been looking for a while and didn't come across any answers yet.


texts: BuildOverviews() not supported for VRT dataset

texts:
I am trying to build overviews in a VRT dataset with python GDAL.

I try this code:

out_vrt = os.path.join(out_dir, 'mosaic.vrt')
ds = gdal.BuildVRT(out_vrt, [raster1.tif, raster2.tif])
factors = [128, 256, 512]
gdal.SetConfigOption('COMPRESS_OVERVIEW', 'DEFLATE')
ds.BuildOverviews("AVERAGE", factors)

, but I get this error:

ERROR 6: /home/umberto/Documents/jrc/FB_Data4Good/mosaic/test/vrt/mosaic.vrt: BuildOverviews() not supported for this dataset.

The VRT mosaic is created, but the overviews are not.

As far as I know, gdaladdo should work with VRTs, but I am not sure if gdal.BuildVRT() is actually using it.

Any idea?


texts: Iterators on the same FeatureCollection provide different results

texts:
I have a WFSDataStore, which is based on Geoserver. The data are coming from one Database table (with columns ID (unique), NAME, YEAR)

I created a FeatureCollection with following code:

FeatureCollection<SimpleFeatureType,SimpleFeature> fc=(FeatureCollection<SimpleFeatureType,SimpleFeature>)WFSFeatureSource.getFeatures(query);

Then I created 3 Iterators on that fc (I wouldn't do it in real code, I just did it to understand the error I got):

FeatureIterator<SimpleFeature> itr1=fc.features();
FeatureIterator<SimpleFeature> itr2=fc.features();
FeatureIterator<SimpleFeature> itr3=fc.features();

When I try to get the attribute values, I get different results depending on the iterator:

while (itr1.hasNext()) {
    SimpleFeature ft1 = (SimpleFeature) itr1.next();
    SimpleFeature ft2 = (SimpleFeature) itr2.next();
    SimpleFeature ft3 = (SimpleFeature) itr3.next();
    System.out.println(" Itr1: " + ft1.getAttribute("ID") + ft1.getAttribute("YEAR"));
    System.out.println(" Itr2: " + ft2.getAttribute("ID") + ft2.getAttribute("YEAR"));
    System.out.println(" Itr3: " + ft3.getAttribute("ID") + ft3.getAttribute("YEAR"));
}

Output:

Itr1: 12458 null
Itr2: 12458 null
Itr3: 12458 2017

Does anyone have an idea, how it can be?

Is that a problem of our geoserver configuration may be?


texts:
I have a shapefile with points and I want to create square buffers around my points 6x6. My coordinate system is in meters. I use R and for the buffers I use:

bufers=gBuffer(file, width = 6,capStyle="SQUARE",byid=TRUE).

When I load the buffers in QGIS, I use the measure line plugin for my windows and instead of 6x6 I see that they are 11X11.

Is that normal?


texts: How to recode bands in which max value occur using ArcGIS Desktop?

texts:
I have 2 raster images; ras1 (an image with max value) and ras2 (an image with 5 bands). Basically. ras1 image extracted from ras2 images.

I'd like to find bands from ras2 in which max value occur and recode them 1 to 5. For example, ras2 has band1, b2, b3, b4, b5. if max value of ras1 occur in band 3, 3 return to output, if max of ras1 occur in band 4, 4 return to output etc.

Could you give me suggestion what tool in Arcgis has ability to do this?

I used inlist tool, but it gave me the values (not the number 1-5).


texts: Installing Anaconda for Python 3.6 without wrecking Python 2.7

texts: I have both Python 2.7 with ArcGIS 10.5 Desktop and Python 3.6 with ArcGIS Pro. A few weeks back, I tried to install Anaconda (w/ Conda, Spyder, Jupyter, etc..) and it somehow overwrote my system settings for Python 2.7 (as well as pip). I noticed that I could no longer use pip or import stuff in 2.7. I uninstalled Anaconda and the updated version of pip (it somehow updated the pip version) and everything in 2.7 worked again. Any suggestions on reinstalling Anaconda without wrecking my Python 2.7 development environment again?


texts: Raster reclassification in QGIS or other open source software with automatically calculated breaks

texts:
I've been searching for a way to automatically calculate breaks to reclassify a raster using QGIS or another open source software, in a similar way it can be done in ArcGIS (Slice tool). I've searched and I found some related questions like QGIS Natural Breaks with raster file, Apply “natural breaks” classification to a raster symbology in QGIS? and calculation of breaks for a quantile classification in GRASS (Raster style based on percentile QGIS). The questions related to natural breaks didn´t have an answer. I've also checked some R packages like classInt, but it seems there are problems for large rasters, as it's shown in this post in the r-sig-geo forum (Jenks classification for raster representation).

I've recently been reading the book QGIS 2 Cookbook (a good book by the way) and they say that the SAGA Cluster analysis for grids algorithm could be used to reclassify a raster. I haven´t found many examples of using this algorithm for raster reclassification, apart from the one in the book Geomorphometry, in the chapter related to Geomorphometry with SAGA and in this presentation (Digital Soil Mapping in SAGA GIS).

What I am trying to do is to convert a quantitative map in a qualitative map, let's say I have a fire risk raster and I want to reclassify it to low, medium, high and very high classes.

Is this algorithm suitable for this kind of raster reclassification?

Are there other ways of automatically calculate the classes breaks using QGIS or another open source software with the Jenks Natural Breaks algorithm or other algorithms different from Equal Interval or Continuous?


texts: Error while performing an intersection of a point and a multipolygon

texts:
I'm willing to build a tiny service that given a point location, it searches within a multipolygon, which polygon object it is intersecting. My current development is defined as this:

from osgeo import ogr, osr, gdal
import geopandas as gpd
import shapely
import json
import requests
import zipfile
import io

'''
Returns the city, state of a certain geolocation in Brazil
'''

debbuging_mode = True

if debbuging_mode:
    ogr.UseExceptions()
    osr.UseExceptions()
    gdal.UseExceptions()

def main(lat=None, lon=None):
    # Variable setup
    if (lat == None) or (lon == None):
        # Using a generic city for demonstration
        with open('static/test/ponta_pora.json') as json_file:
            data = json.load(json_file)
            for key in data:
                lat = data['latDec']
                lon = data['lonDec']
    ibgeCode = city_name = state = None

    # Define spatial reference (WGS84)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)

    # Create ogr point object
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(x=float(lon),y=float(lat), z=float(0.0)) # Forcing datatype is mandatory

    shp_file = 'static/datasets/ibge/2018/UF/BRUFE250GC_SIR.shp'
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataSource = driver.Open(shp_file,0)
    layer = dataSource.GetLayer(0) # Single layer shapefile by definition

    intersection = None
    layer.Intersection(point, intersection)

    print(intersection)
    print('Intersection: '+intersection.ExportToWkt())

    return lat, lon, ibgeCode, city_name, state

The shapefile is a multipolygon with its parts being all cities in Brazil, and the Json provides a valid coordinate in WGS84.

However, I am getting an error when trying to intersect:

Traceback (most recent call last):
  File "/Users/ricardobarroslourenco/PycharmProjects/brazil_city_search/bcs/__init__.py", line 88, in <module>
    main()
  File "/Users/ricardobarroslourenco/PycharmProjects/brazil_city_search/bcs/__init__.py", line 48, in main
    layer.Intersection(point, intersection)
  File "/anaconda3/envs/brazil_city_search/lib/python3.7/site-packages/osgeo/ogr.py", line 2555, in Intersection
    return _ogr.Layer_Intersection(self, *args, **kwargs)
TypeError: in method 'Layer_Intersection', argument 2 of type 'OGRLayerShadow *'

Any hint on what this error might be? Note: The final print is just for debuggig, ideally I would be parsing the intersected plygon object and returning its properties.


texts:
I have an .mxd with a couple dozen layers each pointing to the same feature class and joined table. Each layer is symbolized by attributes in the joined table and have definition queries using the joined table.

In this year's update of the database, our client has requested that we change the name of the joined table - new records have been added and edited, but the schema is unchanged. And since the customer is always right, the sales guys said "of course!"

Now, since the table name has changed, everything is broken. Is there a way to repath all of these joins and definition queries to the new joined table name? "Set Data Source(s)" does not appear to do what I need it to here.


texts: Zooming/panning map in QGIS layout independently from frame?

texts:
I am coming from ArcMap. Working in the layout view, I had buttons to zoom/pan data and other buttons to zoom/pan the view itself. In QGIS, I can't seem to be able to achieve this. The only buttons I have are for zooming/panning the view not the data(the same happens when I use the mouse scroll wheel).

Any way to achieve this?


texts: Configuring setup to allow hosted service to be accessible outside network?

texts:
I have services hosted on my GIS server they are being consumed in maps published to ArcGIS online and also Collector. When i am connected to my network at work, everything works perfect. However when i get off my work network, it all falls apart. I will try to access my services via AGOL or Collector only to get 'services not available. NSURL Error -10003". So somehow my services are not available out side my network

i thought by having the webaptor install on my webserver was supposed to take care of this, and make my services available, but it doesnt seem to be doing that. So what is it about my set up that is off.

Machine 1:
Windows Server 2012R | ArcGIS for server 10.3.1 | Windows Firewall is on and i have the following rules in place: Inbound - Local ports 80, 6080, 4000-4003; Remote Port: All Ports | Outbound - Local ports 80, 6080, 4000-4003; Remote Port: 80, 6080, 4000-4003

Machine 2:
Windows Server 2012R | ArcGIS webadapter 10.3.1 | Windows Firewall is on and i have the following rules in place: Inbound - Local ports 80, 6080, 4000-4003; Remote Port: All Ports | Outbound - Local ports 80, 6080, 4000-4003; Remote Port: 80, 6080, 4000-4003

So what gives why can i not access my hosted features outside my network?


texts: Is it possible to create multiple band using ee.Image.expression?

texts:
I'm using the gee API in a python environment.
I'd like to create several composit to add to my image using the expression function. In the doc they specify :

If the result of an expression is a single band, it can be assigned a name using the '=' operator (e.g.: x = a + b).

So could it be possible to create several bands ?
#just to illustrate
calc = "gfc = (A<={0})*((C==1)*50 + (C==0)*30); " #band 1
    calc += "A+B ;"         #band 2
    calc += "C*B; "        #band 3

    calc = calc.format(threshold)

    bands = {
        'A': clip_dataset.select('treecover2000'),
        'B': clip_dataset.select('lossyear').unmask(0),
        'C': clip_dataset.select('gain'),
    }

    gfc = clip_dataset.expression(calc,bands)


texts: Can Ksp Algorithm(Shortest Path Algorithm) have a constraint that the suggested paths must include a particular node or link

texts:
I am using pgr_ksp algorithm to find out the shortest path between a startnode and endnode by creating a network topology between all the nodes. Pgr_ksp algorithm very well gives me all the shortest paths taking the cost in consideration. But now, I also want to add a constraint that the path its giving between startnode and endnode must include or exclude a particular node or link.   Below code works perfect. But how and where can I add a constraint such that the path will include a node for example 'nodeA'?

SELECT d.*, f.linkname ,SUM(cost) OVER (PARTITION BY path_id) AS aggregation
FROM pgr_ksp('Select id, source, target, airmiles as cost from f',
source,
target,
 2) AS d
INNER JOIN f ON d.edge= f.id
ORDER BY aggregation, path_id, path_seq;


texts: How to get a QGIS useable style file from OS .xml and .xsl files?

texts: Ordnance survey provide some style files to support their vectormaplocal product. The style files appear to come in pairs, similarly named but one with a .xml extension and the other .xls   The two files together seem to be designed to define the style but my problem is how do I get the two into a format that QGIS can handle?


texts: Copying attributes from one polygon layer to another?

texts:
I have a problem that I can't seem to get my head around. I have two polygon layers:

Polygon A - is a subset of polygon B with the same fields and has identical polygons to Polygon B
Polygon B - has the attribute data I want to be in Polygon A

How can this be done? I tried the QGIS tool "Join Attributes by Location" but as some of the polygons are within others, it tends to link to the first intersect it finds (the outer polygon).


texts: TypeError: unhashable type: 'dict' when using shapefile with arcpy_metadata module

texts:
I am having an issue with the arcpy_metadata module. I have a simple shapefile I am trying to load into the MetadataEditor:

import arcpy, os, sys,env
import arcpy_metadata as md
myfeatures = "C:/Workspace/Sandbox/MetadataProject/data/model_domains"
arcpy.env.workspace = myfeatures
list = arcpy.ListFeatureClasses()
for item in list:
    metadata = md.MetadataEditor(item)

When I run this, I get the following error:

TypeError: unhashable type: 'dict'

How can I resolve this?

It also seems that since I am running 64-bit version, that features from a .gdb won't work becuase they require XSLTransform_conversion, which is only supported in 32-bit.


texts: Use a non-GDB SDO_GEOMETRY table in a map service?

texts:
I have a table with an SDO_GEOMETRY column (points) in an external system.

The external system is an Oracle 12c database; it is not a geodatabase.

I want to serve up the SDO_GEOMETRY table to a map service/ArcGIS Server 10.7.1.

I want to display & query the layer in a web map which is based on a map service. I do not need to edit features.

Is it possible to do this, even though the table is not registered with a geodatabase?


texts:
I have a raster map that needs to be defined in 2 different coordinate systems.
I added control points to my unprojected map and rectified it. It's OK for the first coordinate system. (European datum ed50 30 3).

Now, I want to project this raster to ed50 30 6. So right click to my raster under catalog | properties | spatial reference tools, and I choose new projected system and click ok.

When I restart ArcGIS, adding my new raster, the coordinate system doesn't change. Projection is listed as (ed50 30 6) but the map location doesn't change (I can understand from grids on the map it still has same coordinates)

Am I skipping some steps?


texts:
I have a PHP endpoint/service that receives Lat/Long data from a 3rd party. The data is then processed in PHP code and saved into a MySQL database. This process works fine without any issues.

I need to run some spatial queries against the data, like does the source point fall inside a polygon (polygon source is GeoJSON), what is the distance from the point to a line (line source is GeoJSON) and distance to the end of a line (following the line segments).

I found TurfJS can cover all these queries, BUT TurfJS is obviously client side JavaScript functions, so I'm unable to call the Turf functions from PHP, retrieve the output values to save it in MySQL.

Any suggestions?


texts:
I recently used the GPSlogger app and created KML files from my routes. I haven't been able to successfully convert them to shape files.

Does anyone know a good downloadable program that converts file types?


texts:
I am new to GIS.
I am using Mapinfo V17.0.2
I am creating a cemetery plot. I have created many rows of polygons to shape the desired areas for burials.
My table looks a little like this;

ID | Lot | Sector | Long | Lat |Cemetery

1  | A1 | JE | 664915.21 | 6334263.07 | GA

I have finished manually entering Row and Lot into one of my table columns [Lot]; ie. A1, A2....W23,W24. ( If there is an easier way I'd be much appreciative)
Using a raster image of the area I Extracted the Long and Lat.

Now I needed different 'Sector's within, Like Anglican, Jewish etc..
'AN' & 'JE' respectively.
Now I entered the JE manually, like 10 entries. I knew I could export these to excel, .csv format and then autofill the remainder with 'AN'... I went to Data > Filter and showed the top entries with JE so i could autofill the rest with AN. Now what I didn't know was by doing this when I opened the table back into Mapinfo, all the layer objects looked correct but my data was all over the place. where A1 should of been was now J26 , this was because I filtered the data in excel.

How can I resolve this, or do I need to re-start my manual Row and Lot entries again?

Looking at my structure and how I started to create my polygons, I am not sure how to fix this. As I created polygons for the cemeteries that are existing and then creating polygons for empty plots so my table strucure already wasn't in order of A1..A2..A3 etc. It would of looked more like; J3..J4..E5..E12..A1
Another day just manually inputting the Row and Lot back in I'd say. :/


texts: I signed into GEE code editor today to find that all of the scripts I've written and saved are now gone. It just says "No accessible repositories" under owner/writer/reader drop down panel. All of my assets are still there though. What should I do?


texts: latitude longitude search in ALL state or county gis sites

texts:
I'm working in geography field now and
using local gis sites for different states and counties.

However, they all seem different, and it takes time to
learn each new one as required.

For example, now I am using the Maine gis at
http://mapserver.maine.gov/basemap/index.html.
I found the way to search by lat/long on the
last county gis I was on in colorado, after some
searching for it.  However, I cannot find such a
feature on the maine one yet, and it's starting to take
too much time.

Can anyone help me find the lat/long coords search function
on that site, or else refer me to a more global gis search site
that will provide property parcels, addresses, lat/long,
so that I won't have to learn each state or county gis site
individually?

I am new to gis so I am a bit naive about the whole thing.


texts: EU DEM in GRS80 EVRS2000 with geoid EGG08 --> Make it comparable to Black Sea 1975 for Romania

texts:
I got contour lines for western Romania in vertical system Black Sea 1975 which is based on elevation data with 100 m resolution. Now I would like to compare that to the EU DEM with higher resolution. How can I make those two comparable? Is there a certain Z value that I can use to lower or elevate the EU DEM?

I want to compare two DEMs with two different geodetic systems or "reference sea levels". And those two reference sea levels are "Black Sea 1975" and "GRS 80 EVRS2000 with geoid EGG08". And my question is how I can make those DEMs comparable since they - as I assume - have a different 0. The 0 of "Black Sea 1975" is somewhere in the Black sea and I don't know where the 0 of GRS80 is.

Here is a famous case about using different reference sea levels. I hope that makes my question very clear now..


texts: Reduce frequency of vector points based on attribute or spacing

texts:
I have a point shapefile with 5000+ points.  The points represent centroids of a flood model mesh output and I want to display the flow direction attribute in the form of an "arrow" point style.  I want to display for purposes of printed maps.

The points are too dense to allow practical presentation.  I would have preference for retaining all data, but ultimately I need a way to "thin out" the points so that I only retain points that are outside a proximity to another point and/or is significantly different in the direction attribute to the adjacent points.

How can I achieve this / can I achieve this in QGIS?


texts:
I'm using ArcMap/ArcScene 10.4.1

I am trying to extract z-values out of the file in the picture below, but can't find any way to do it. The height is different for the points, so it has z-values, but how do I extract them?


texts: OpenLayers.Filter.Comparison value charecter encoding

texts:
I am using vector layer with WFS protocol

var adressVector = new OpenLayers.Layer.Vector("Adress", {
  strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: new OpenLayers.Protocol.WFS({
    url: "http://localhost:8080/geoserver/wfs",
    featurePrefix:"compa",
    featureType: "adress_location",
    featureNS: "http://localhost:8080/location",
    srsName: "EPSG:4326",
    geometryName: "the_geom",
    version: "1.0.0"
  })
});

I want to get addresses with some filter and adding filters property fo vector layer:

var adressVector = new OpenLayers.Layer.Vector("Adress", {
  strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: new OpenLayers.Protocol.WFS({
    url: "http://localhost:8080/geoserver/wfs",
    featurePrefix:"compa",
    featureType: "adress_location",
    featureNS: "http://localhost:8080/location",
    srsName: "EPSG:4326",
    geometryName: "the_geom",
    version: "1.0.0",
    "filter": new OpenLayers.Filter.Logical({
          type: OpenLayers.Filter.Logical.OR,
          filters: [
                    new OpenLayers.Filter.Comparison({
                        "type": OpenLayers.Filter.Comparison.EQUAL_TO,
                        "property": 'name',
                        "value": "BÖRT De Sem"
                    })
                ]
          })
  })
});

But this filter creates a WFS XML like this.

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<wfs:Query typeName="compa:adress_location" srsName="EPSG:4326">
  <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <ogc:Or>
       <ogc:PropertyIsEqualTo matchCase="true">
          <ogc:PropertyName>name</ogc:PropertyName>
          <ogc:Literal>B�RT De Sem</ogc:Literal>
       </ogc:PropertyIsEqualTo>
   </ogc:Or>
</ogc:Filter>
</wfs:Query>
</wfs:GetFeature>

I mean my vaue property is ugly charectes.

"BÖRT De Sem" is going "B�RT De Sem"

So no feature found.

How can I solve this problem?


texts: Finding length and width of polygon using QGIS

texts:
How can I find polygon's length and width in QGIS (regular polygon and irregular polygon)?

I have area and perimeter but not length and width.


texts: I want to restrict where the user can pan on the map, but I can't seem to find any resources that tell how to set the max extent of a map view in OpenLayers 3. I know there was a solution in OpenLayers 2. Is it possible in the new version?


texts: How to split a line at a vertex in QGIS 2.2?

texts: Is there a plugin out there that will split a polyline at vertexes?  My arcgis license isn't high enough to run the tool.  I'm using 2.2.


texts: Checking if road passes through a trapezoid?

texts:
I have a straight road defined by 2 points with latitude and longitude.
I also have a trapezoid defined by 4 points with latitude and longitude.
How do I check if the road passes through that trapezoid?