texts: Geographically Weighted Regression fails QGIS

texts:
So I'm trying to run GWR on a dataset in QGIS.

It starts and then I get these messages at the end:

"Converting outputs
'NoneType' object has no attribute 'dataProvider' See log for more details"

Other error messages in the log that look relevant are:

GDAL: Export Raster: could not initialize data objects

GDAL: Export Raster: could not initialize data objects
Error: Grid file could not be opened.
Error: executing tool [GDAL: Export Raster]

Other possibly relevant details:

Both of the variables are Integers

The dataset has ~9800 entries

Any ideas about what to try or ideas about why its failing?


texts: How to reproduce this type of server-side authentication?

texts: A Server Connection that I have is username and password protected, and so this prompt pops up when you open the mxd. I am trying to understand the system we have here better, and would like to be pointed towards a resource that explains how to set this up. I am unsure if this is an ArcGIS or a server feature.


texts: Flood fill formula for Voronoi cells on QGIS?

texts:
At the moment, I'm developing a map of a fictional West Indian nation with QGIS. Everything's gone as well as can be expected with this project--except for a shortcoming that stumped me on Thursday afternoon.

The latest part of my project involves grouping areas on a Voronoi chart, based on attribute. For instance, if I choose at least two cells and set their value as 1, then every cell touching their sides must have 2 as their value; those touching the "2" cells receive 3; and so forth until they reach the chart bounds. This method, known as a "flood fill", is used in procedural generation; see it in action in this Imgur album.

When trying to select empty touching cells with this expression:

CASE
WHEN "Dist Lvl" = 1
THEN geomtouches(buffer(centroid($geometry), 100))
END

no results appear before I can apply "2" to those.

Before I continue, I'd like to obtain a better code, or a Python/PyQGIS script, that can do the job faster; I struggled to draft one on paper a while ago. (Although I can efficiently use QGIS otherwise, I sometimes have trouble with the PyQGIS bit, a couple of examples aside.)

At least reading so many of these posts here in the past few years--this here is my first one--have really helped me on my way as a geofiction cartographer. So now it's your turn to help me out.


texts: Creating condition for mosaicing images in Google Earth Engine JavaScript API

texts:
I am trying to create a mosaic of a several classified images using as criteria that every pixel has to be the same class at least for three images.
A similar process can be done using the .min() or .max() functions, that will return the pixel corresponding either to the min or max value. Example:
// Collection of classified images
var collection = ee.ImageCollection([image1, image2, image3, image4])

// Function to select a specific class
var selectClass = function(img) {
    var select = img.eq(2); //Select class 2
    var maskImg = img.mask(select); //Mask other classes and let only class 2
    return maskImg.unmask(0) //Assign a no data value to masked pixels
}

// Create new collection with only class 2
var collectionNew = collection.map(selectClass)

// Create mosaic using .min(). If there is at least one stacked pixel with the value 0,
// then it will return this instead of class 2.
var mosaic = collectionNew.min()

Now, what I want is to create a mosaic with the condition of returning the pixel value of class 2 if there is three or more pixels with the same value, otherwise return 0.
Is this possible in GEE using the JavaScript API?


texts: Unable to connect to PostgreSQL database using schemas

texts:
I am using QGIS 2.0 with PostgreSQL 10 and PostGIS 2.4.

I have multiple users who need to access their vector data from the database. All the users have the same vector layers on their project, but need to create their own features and store them individually on the database, with the user being able to access only his/ her own tables.

I created a database user_markings and imported all the relevant vector layers into it using DB Manager. I cloned the public schema of this database using a script into multiple schemas, one for each user. I granted all privileges for a user to his/ her schema, and assigned the user login a path to the schema.

When I attempt to connect to the database as one of the users, I get the following error: Erronous query: SELECT postgis_version() returned 7[ERROR: function postgis_version() does not exist].

I am able to login as a Superuser. I have double checked the postgis and postgis_topology extensions in my database. Also as I read in one of the posts, I have assigned a path to postgis in the database properties.


texts: I am new to GIS development and was wondering if there are tutorials available for development in ArcGIS 10. I have some experience writing simple scripts in python and fooling around in model builder, and I want to take the dive and make an effort to become a real programmer.  I need some more info to get a better grasp on doing some real GIS programming to make some custom tools in ArcGIS 10. Where/how did you guys learn all your GIS development skills? Any book/web site recommendations? Just an FYI I will being doing my coding in visual studio using the .net framework.


texts: Extracting time series Soil moisture values using point data in Google Earth Engine?

texts:
I am writing a script to extract the time series soil moisture values from SMAP soil moisture data using point shapefile and plot it on a line graph in the google earth engine. As SMAP having the daily data coverage but I am getting the pixel values three days interval. Please can anyone tell me where is the wrong with my code and addition to that I want to extract the time information also along with date and export it as csv to drive,
Below code, i am working
#Point shape file i uploaded, if i can add manually coordinate, please tell me (Not from geomarty addition)
var shp=ee.FeatureCollection('users/SWAT/C01');

// which dataset + where + when
var imageCollection = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture')
                .filterBounds(shp)
                .filterDate('2018-01-01','2020-08-30');

// select bands -> ssm = surface soil moisture (%)
var soilmoisture = imageCollection.select('smp');

// gget scale
var scale = imageCollection.first().projection().nominalScale();

// Take average
var imageCollection = imageCollection.map(function(image){
  return ee.Image(image.setMulti(image.reduceRegion(ee.Reducer.mean(), shp, scale)));
});

// Get the mean from the dictionary and print it.
var means = imageCollection.aggregate_array('smp');
print('Soil Moisture (%)', means);
// put in a chart
print(ui.Chart.image.seriesByRegion(imageCollection, shp, ee.Reducer.mean(), 'smp', scale));

// Get the mean from the dictionary and print it.
var meanDict = soilmoisture.mean().reduceRegion(ee.Reducer.mean(), shp, scale);
print('Soil Moisture (%)', meanDict);

Map.centerObject(shp, 10);
Map.addLayer(soilmoisture, {min: 0, max: 1}, 'soilmoister imageCollection');


texts:
I've previously converted simple shapefiles to WKT so that I can store them in RDMS that has spatial capabilities. I'm a new user to QGIS, which I'm using to browse a shapefile. But now I have a shapefile like this:

Not only are there polygons, but there are labels on those polygons. Really what I'd like is to strip off the polygons and save them with their associated label. So something like:

row1 ("POLYGON((1 1,5 1,5 5,1 5,1 1))) 'high'
row2 ("POLYGON((2 2,5 1,5 5,1 5,2 2))) 'medium'
row3 ("POLYGON((5 5,5 1,7 7,1 5,5 5))) 'high'
row3 ("POLYGON((8 8,5 1,5 5,1 7,8 8))) 'low'

I'm a little stumped at the best way to creatively slice up my shapefile, ideally by some of the attributes defined in the accompanying dbf file. Ideally I'd like to:

export, as WKT, all of the polygons with 'high'
export, as WKT, all of the polygons with 'low'

etc, etc.


texts: How to install or build GDAL 3+ on Debian buster?

texts:
I'd like to install GDAL 3+ on Debian stable. Is this not possible? I've read GDAL docs and looked at Debian package tracker.
With the default source.list for Debian buster (current stable as of October 2020), only GDAL 2.4 is available. Any solution recommending adding the unstable source in this won't work as you'd have to upgrade to Debian unstable.
EDIT: conda is a solution but doesn't work for me. This isn't the best solution because someone might not want to use a conda-based environment.


texts: Making Hierarchy "less strong" in ArcGIS Network Analyst?

texts:
I want to know a way where the Hierarchy in my network analyst doesn't influence that much.

What happens is that when making shortest route it makes a 25-mile route, with the hierarchy it chooses to make a 40-mile route, choosing a principal road that goes a lot further than I need.

How can I make my Hierarchy more customizable, so it picks the better roads but without making so larger routes?


texts:
i tried out this recipe from the GDAL/OGR cookbook:

from osgeo import gdal, ogr
# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = -9999
# Filename of input OGR file
vector_fn = 'test.shp'
# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'
# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res,        1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)
# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])

I recieve a rectangular raster in the area were my polygon is located. I want to mask the raster exactly to the boundaries of a multipolygon shapefile. Is there an option to GetGeometry from the shapefile and use this as parameter instead of(x, y).


texts: Drawing line that connects pins in ArcGIS Online Story Map?

texts:
I need to draw a line in the streets connecting the pins. I have a public account.

Does any one know how to do it?

Attached is what I'm trying to accomplish.


texts: Complete list of UK postcodes plus their city, county, uk country

texts: Is there a free data set of all current UK postcodes (about 1.7 million-ish) plus their city name, county and UK country (e.g. Liverpool, Merseyside, England) and ideally lat, long, east, northings?


texts: How to best remove rgb values outside of aerial image?

texts:
I have just converted RGB-TIF images to RGB-ECW for reasons of space and performance in ArcGIS. However, after the conversion, there are relics outside of the images which are not nodata values but have the values 0,1,0 or 1,1,0 or 1,2,1 etc. All very low values which are shown as black remnants after I set nodata values to transparent. (See attatchment)

I have ArcInfo and FME at hand and would like to either change these rgb values to 0,0,0 or clip all values outside of the aerial extents. Could anyone offer a tip on the most efficient way to achieve this?

Image remnants:


texts: Control legends from multiple data frames simulataneously

texts:
Is it possible to control the symbology of all legends contained in multiple data-frames (but on a single layout) simultaneously?

I want them to all have the same style, font etc.


texts: GeoServer - Sometimes bounding box property is missing following CQL Filter - why?

texts:
Usually when I run a WFS GetFeature query on GeoServer and request JSON the results conform to this format:

{
"type":"FeatureCollection",
"features": [..features here...],
"crs":{"type":"EPSG",
    "properties":{"code":"4326"}
      },
    "bbox":[0,52,0,52]

}

Yet I've started to find that when I run certain CQL filters on the GetFeature query, the bbox property is missing in the results.

For example:

madeupserver/geoserver/layergroup/ows?service=WFS&version=1.0.0&maxFeatures=5000&request=GetFeature&typeName=dbossh:fullflat&outputFormat=json&CQL_FILTER=name like '%25a%25' and owner like '%25bob smith%25'

returns the full JSON result set with the bbox, but

madeupserver/geoserver/layergroup/ows?service=WFS&version=1.0.0&maxFeatures=5000&request=GetFeature&typeName=dbossh:fullflat&outputFormat=json&CQL_FILTER=name like '%25ab%25' and owner like '%25bob smith%25'

(note the addition of a single letter to the name like part of the query).

returns a result set which is missing the bbox.

I'm really curious now as to why the bbox would be missing. Can anyone help me understand what the difference is that would cause GeoServer to leave out the bbox?


texts: R st_centroid GEOS error Unknown WKB type 12

texts:
I’m trying to calculate centroids of polygons of postcode areas (wfs shapefile), but I’m getting the following error message:
Error in CPL_geos_op(“centroid”, x, numeric(0), integer(0), numeric(0), : Evaluation error: ParseException: Unknown WKB type 12.

Is there a way to fix this? I have no idea what this means. I assume there is something unusual with the formatting of the geometry column(?).
> sessionInfo()
R version 4.0.4 (2021-02-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] ows4R_0.1-5     geometa_0.6-3   forcats_0.5.1   stringr_1.4.0   dplyr_1.0.5
 [6] purrr_0.3.4     readr_1.4.0     tidyr_1.1.3     tibble_3.1.0    ggplot2_3.3.3
[11] tidyverse_1.3.0 httr_1.4.2      sf_0.9-7        sp_1.4-5        RPostgres_1.3.1

Reproducible example (get data from WFS server)
library(tidyverse)
library(sf)
library(httr)
library(ows4R)

# get postcode data ####

# following info from
# https://inbo.github.io/tutorials/tutorials/spatial_wfs_services/

# connect to Statistics Finland geoserver ####
wfs_regions <- "http://geo.stat.fi/geoserver/postialue/wfs"
regions_client <- WFSClient$new(wfs_regions,
                                serviceVersion = "2.0.0")
regions_client$getFeatureTypes(pretty = TRUE)

url <- parse_url(wfs_regions)
url$query <- list(service = "wfs",
                  request = "GetFeature",
                  typename = "postialue:pno_2018",
                  srsName = "EPSG:3067"
                  )

request <- build_url(url)

# grab data ####
postal_codes <- read_sf(request)
# st_geometry(postal_codes)

ggplot(postal_codes) +
  geom_sf()

# subset to smaller area, e.g. ####
pc.helsinki <- postal_codes[postal_codes$kunta == "091",]

ggplot(pc.helsinki) +
  geom_sf()

# centroids of postcode polygons ####
pc.helsinki$centroid <- st_centroid(pc.helsinki$geom)


texts: Extracting convex_hull polygon from CSVwith scipy.spatial ConvexHull and output it to a shapefile

texts:
I'm looking for a Python solution for extract a convex_hull polygon from quite big CSV (nearly 150 millions points) and return back the polygon to a shapefile.
After reading the CSV with pandas (pd.read_csv), I give a change to  scipy.spatial Convex_Hull due to very quick response time. And effectively I get the convex_hull envelope correctly estimates.
But I'm now in the hurdles, because I've not a solution to  export the result (a scipy.spatial.qhull.ConvexHull ). If I understand correctly, I need to extract from the hull structure the points or vertices, because no polygon is available.
Any hints?


texts:
I am querying a list of the worlds parks with something like this query:

http://api.geonames.org/search?featureCode=PRK&maxRows=10&username=demo

but is it possible to query for the type of park it is (county/national/state/city) or also to break it down by state and country, or just query within one country?

Is that asking too much or is that possible to do with that API?
How could I get those values?  And how do I get the total number of parks they have?


texts: Skipping outputs from the QGIS Batch Processing?

texts:
I am doing a lot of processing in batch with the r.sun model from GRASS 7.2.1 in QGIS 2.18.9.
I am using the r.sun model with intention to reuse for further raster calculations the Irradiation raster and the Diffuse irradiation raster. The model computes also the Ground reflected irradiation, the Global irradiance and the Insolation time rasters, which I don't need for my further analysis.
So, I need just 2 of 5 output calculations for making further analysis.

I am interested in if I can skip somehow this outputs. I know I can put them easily all in one folder/TMP with autofill and then delete them after, but I would love just to skip their creation by default to spare the processing time as it is really long (30min per computation).

I can manage some simple code lines and the only idea I got is to export this batch process in JSON and set the skipping of this data creation there and relaunch it from the console.
Any suggestions about how to do it or some better idea for someone who is not able to program it by scratch?


texts: Classifier with SVM algorithm for sentinel 1A image

texts:
I am classifying rice for the Red River Delta with 2 class:
class 0: Rice
class 1: Not rice

I used svm classifier algorithm. However, after running the code, the classifier operated for a long time without results.
Here is my code :
Link:
https://code.earthengine.google.com/f70f6ae7d54a428829a84cbc9d27ea3b

Code:

//Load 1 shapefile red river delta
var dbsh_shape = ee.FeatureCollection('users/chungpd/DBSH_Shape');
Map.centerObject(dbsh_shape);
Map.addLayer(dbsh_shape);

// Use these bands for prediction.
var bands = ['VV', 'VH'];

// Load an imagecollection over a portion of Red river delta
var img1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterDate('2017-01-03', '2017-01-18')
  .filterBounds(dbsh_shape);

// Cut img1 according to "dbsh_shape" shapefile
var img2 = img1.map(function(image) {return image.clip(dbsh_shape);});

// Mosaic img2
var image = img2.mosaic();
print (image);

// Load training polygons from a Fusion Table.
// The 'class' property stores known class labels.
var polygons = ee.FeatureCollection('ft:1fxOJdG56CUWXy7MWi4KPeCAvKYC4YHFzMbRPU7S3');

// Get the values for all pixels in each polygon in the training.
var training = image.sampleRegions({
  // Get the sample from the polygons FeatureCollection.
  collection: polygons,
  // Keep this list of properties from the polygons.
  properties: ['class'],
  // Set the scale to get Landsat pixels in the polygons.
  scale: 10
});

// Create an SVM classifier with custom parameters.
var classifier = ee.Classifier.svm({
  kernelType: 'RBF',
  gamma: 0.5,
  cost: 10
});

// Train the classifier.
var trained = classifier.train(training, 'class', bands);

// Classify the image.
var classified = image.classify(trained);

// Create a palette to display the classes.
var palette =['00FF00', 'DD0000'];

// Display the classification result and the input image.
//Map.addLayer(image, {bands: ['VV', 'VH'], max: 0.5, gamma: 2});
//Map.addLayer(classified, {min: 0, max: 10, palette: palette}, 'Isrice');
Map.setCenter(105.96999, 20.90339, 7);
//Map.centerObject(image, 10);
Map.addLayer(image, {bands: ['VV', 'VH'], max: 0.4}, 'image');
Map.addLayer(classified, {min: 0, max: 1, palette: ['00FF00', 'DD0000']},
  'classification');

This is the error:

classification error: Layer error: Computation timed out.


texts:
My friend and I, we are developing a plugin for qgis but we have a drawback when doing the multibuffer we do not understand the methods well
and the code we are using is the following

#Generación de los buffer sobre la capa vías
        # Se leen los vectores que esten el la lista de capas cargada en QGis
        layers = self.iface.legendInterface().layers()

    # Selecciona la capa vial para hacer multibuffer (capa de vías)
    selectedLayerIndex = self.dlg.capavias.currentIndex()
    selectedLayer = layers[selectedLayerIndex]

    # Generar el buffer de la capa vias
    vias = selectedLayer

    QgsGeometryAnalyzer().buffer(vias, "C:/Users/toshiba/Downloads/buffer_vias100.shp", 100, False, True, -1)
    QgsGeometryAnalyzer().buffer(vias, "C:/Users/toshiba/Downloads/buffer_vias200.shp", 200, False, True, -1)
    QgsGeometryAnalyzer().buffer(vias, "C:/Users/toshiba/Downloads/buffer_vias300.shp", 300, False, True, -1)
    QgsGeometryAnalyzer().buffer(vias, "C:/Users/toshiba/Downloads/buffer_vias400.shp", 400, False, True, -1)
    QgsGeometryAnalyzer().buffer(vias, "C:/Users/toshiba/Downloads/buffer_vias500.shp", 500, False, True, -1)
    buffer_vias500 = iface.addVectorLayer("C:/Users/toshiba/Downloads/buffer_vias500.shp", "", "ogr")
    buffer_vias400 = iface.addVectorLayer("C:/Users/toshiba/Downloads/buffer_vias400.shp", "", "ogr")
    buffer_vias300 = iface.addVectorLayer("C:/Users/toshiba/Downloads/buffer_vias300.shp", "", "ogr")
    buffer_vias200 = iface.addVectorLayer("C:/Users/toshiba/Downloads/buffer_vias200.shp", "", "ogr")
    buffer_vias100 = iface.addVectorLayer("C:/Users/toshiba/Downloads/buffer_vias100.shp", "", "ogr")


texts:
I want to create a decimal degree (dd) field from degree minute seconds (DMS) as a script tool.

The CSV looks like this.

My code looks like this:

import arcpy
inTable = arcpy.GetParameterAsText(0)
x_deg = float(arcpy.GetParameterAsText(1))
x_min = float(arcpy.GetParameterAsText(2))
x_sec = float(arcpy.GetParameterAsText(3))
y_deg = float(arcpy.GetParameterAsText(4))
y_min = float(arcpy.GetParameterAsText(5))
y_sec = float(arcpy.GetParameterAsText(6))
DD_X = (x_sec/3600) + (x_min/60)+ x_deg
DD_Y = (y_sec/3600) + (y_min/60)+ y_deg
arcpy.AddField_management(inTable, DD_X, "DOUBLE", 20, 10)
arcpy.AddField_management(inTable, DD_Y, "DOUBLE", 20, 10)

I keep getting an error message:

  ValueError: could not convert string to float: 'x_degree'

I feel like this should be straightforward, but perhaps I'm over simplifying.

Can someone point out what I'm doing wrong?


texts: Determining ESRI File Geodatabase release version?

texts:
Is there any method in Python, ArcObjects, or the File GDB API for determining the release version (9.3, 10.0, etc.) of a file geodatabase? There is a need to find all of our shared 10.0 version file geodatabases and convert them back to 9.3 since not everyone as at 10.0 yet.

I know that you can create 9.3 version file geodatabases in 10.0 so that part should be easy enough.


texts: Identifying unknown map projection for image

texts:
I intend to georeference and create a shapefile for some of the features but I'm unable to determine the correct map projection. (Working in QGIS ideally.)

The image is from a scientific study but the map projection is unknown. Does anyone have an insight about this map projection?


texts: Removing n first digits on a field attribute table using QGIS string functions

texts:
I'm looking for a QGIS function to remove the first n digits on a field attribute table using QGIS string function.
For exemple I have the following string in a QGIS attribute table:

  742-004-SRV Superficie : 12780

My aim is to remove the first 12 characters. How to achieve this with QGIS?


texts:
We currently use Pitney Bowes Confirm Connect for our highways inspections. We are trying to use our own WMS feed instead of the default Bing Maps to allow us to have access to information such as adopted highway records.

Following the guidance document you have to write a javascript settings file and drop it into the Confirm folder on our Android device (Samsung Tab 3). The application picks up the copyright text but displays a blank basemap. Can anyone see anything wrong with the JS below?

{
"mapSettings":{
"roads":{
"type": "WMS",
"projection": "EPSG:27700",
"projectionDef": "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs +nadgrids=ostn02",
"attribution": "(c) Crown Copyright",
"units": "m",
"maxExtent": [250000,150000,400000,180000],
"url": "http://council.company.com/getows.ashx",
"parameters": {
"layers": ["highways"],
"url": "mapsource=Banes/AllMaps"
}
}
}
}


texts: Picture Labels with leader line from point without annotation?

texts:
The closest thing I have is using a line call out in the text background, then placing the picture in the line call out border. But the picture tiles throughout. Is what I am trying to do even possible without adding graphics or annotation?

Basically I want to do what Google earth does automatically with the placement of their labels?

Any extensions besides Maplex? Any tools or solutions?


texts: Comparison of Open Source Desktop GIS Packages

texts:
I'm curious as to what comparisons are out there between different Open Source Desktop GIS software packages. I'm aware that there has been a QGIS and gvSIG comparison on this site but I am looking for comparisons between any OpenSource Desktop GIS.

A substantial list of Open Source GIS software can be found in response to the question: What are some Free and Open Source GIS Desktop packages

I don't expect a single comparison with everything in it, but comparative subsets. Also any comprehensive benchmark/testing that may have been performed between any grouping of these products.


texts: ASCII grid (.asc) to reprojected geoTIF in EPSG:3857

texts:
I have an ASCII gridded file with a .prj file along with it. My only goal is to convert the file to a GeoTIFF, reprojected to EPSG:3857.
The prj file only includes the following :
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]]

I found a resource online (https://developers.arcgis.com/javascript/3/jshelp/gcs.html) that matches this prj file with EPSG:4269. Am I understanding this correctly? If so, how do I convert the .asc file to 3857?
I have tried many variations of code but I am no expert. I first tried to convert to TIF :
gdal_translate -a_srs EPSG:4269 in.asc converted.tif -of GTiff

And then reproject to 3857 :
gdalwarp -s_srs EPSG:4269 -t_srs EPSG:3857 converted.tif converted-projected.tif

This results in a projection that is incorrect, but seems close.
What am I doing wrong?


texts: Aggregating/Averaging NDVI across pixels. Arithmetic mean seems bad choice

texts:
I need to compute an aggregate index for vegetation density in certain (large) shapes, to assess desertification. I have NDVI data at a high spatial resolution (~100 NDVI-pixels per aggregation shape, based on MODIS NDVI version 6).

Initially my idea was to compute the arithmetic mean across all pixels in each shape, but that seems to be a bad idea: NDVI is a (unitless) ratio, ranging from -1 to 1, taking approximately the following values:

Clouds/Snow:    ~ <0
Water:          ~-0.1 to 0.1
Soil:           ~ 0.1 to 0.2
Vegetation:     ~ 0.3 to 0.8.
Built-up areas: ~ ?

That scale is not even approximately linear in biomass, area covered in green, or similar natural scales. Taking the arithmetic mean across all pixels in a shape thus does not do what I want, and the resulting aggregation is difficult to intepret. In other words: The arithmetic mean of all pixel-NDVIs in a shape is not the NDVI of that shape. E.g. A change from bare soil (0.1) to concrete (-0.5) in one pixel will affect the average as much as several pixels going from loose to dense vegetation.

Are there any suggestions on how to proceed? In particular I'd like to know what the standard procedure in the remote sensing literature is. Are people actually using arithmetic means, if so, how do they argue? Are there alternatives that are commonly used?


texts: Geoprocessing service to import CAD file

texts:
I need to publish a geoprocessing service to import CAD file in a Web AppBuilder application and add it as a map service.

I tried to use ModelBuilder using the CAD to Geodatabase tool
Made the input cad as model parameter
Passed the output to the CAD to Geodatabase tool

But I can't pass the copy output to the CAD to Geodatabase tool
otherwise the input file will be a static file not from a user upload

I got this methodology as referenced in this GeoNet post.

I tried the tool the person attached but it has its own errors and couldn't resolve due to my lack of Python knowledge.


texts: Adding 2 hours in ArcMap's field calculator using Python

texts:
I am trying to add 2 hours from a particular time, but I do not want to end up with 13:00. I isolated the hour in its own field so the minutes would not be a problem. I am not sure what I am missing. I have switched up the syntax by adding and removing the colons and also capitalizing "if" and "else."

Essentially I just don't want 11am become 13pm instead of 1pm.

I am using the calculate field tool in ArcGIS for Desktop 10.5.1.

if (int(!LT_Hour!) + 2) >12
   x = 1
else
   x = (int(!LT_Hour!)+2)
end if

Any ideas?


texts:
I've made a python function that returns the user logged in a postgis layer:

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def getUser(lay, feature, parent):
    layer=None
    for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyr.name() == lay:
            layer = lyr
            break
    connInfo = QgsDataSourceURI(layer.dataProvider().dataSourceUri()).connectionInfo()
    success, username, password = QgsCredentials.instance().get(connInfo, None, None)
    if success:
        return username

The problem is that the result of the function seems to be ok in the output preview when creating a new calculated field:

But when I check the attributes of any feature of the layer, the field is always NULL:

Not only that but also keeps asking me to insert the credentials for the layer again and again...

The credentials are not stored in the layer so every time I open the map I am prompted to insert them. I did it just once in the beginning until now.

What am I doing wrong?


texts:
I have a problem with the way that QGIS handles transformation parameters between coordinate systems based on the Datum NAD27 and those based on the Datum WGS84. In version 3.12.1 I have found that when loading a layer in NAD27 it asks for the transformation to be used for reprojection. For example, I have data from the mid-1990s that is based on NAD27 which is close to Mexico City. Data is EPSG 26714, UTM NAD27 14N.
Of all the options available, most do not apply because they are designed for Canada or the United States or use NTV grids. I have two options for Mexico, which are:

Inverse of UTM zone 14N +NAD27 to WGS84 (18) World -N Hemisphere -102°W to 96°, Mexico Onshore

Inverse of UTM zone 14N +NAD27 to WGS84 (83) World -N Hemisphere -102°W to 96°, Mexico Offshore GoM -Tampico Area

I choose transformation 18, which has parameters documented by NIMA and is also used by ArcMap, however , I always get this warning message:

I am aware of this question, Reprojecting between NAD27 and WGS 84?

So I wonder What I have to do? Do I rely on Geographic transformation 18? Do I define it manually?


texts:
I'm looking for a global geoapi for the places.

I see Facebook places, google places, simplegeo, factual etc...

Why not create an API that includes all places? Free and easily accessible as OpenStreetMap for the webmapping ?

If people are interested, write me!

Thanks


texts: Decimal Numbers in the Graphical Modeler

texts: How do I set up an input field that accepts decimal numbers in the QGIS graphical modeler?


texts: Removing "noise" from edges of sid Imagery File using ArcGIS Desktop

texts:
I have several .sid imagery files of the state of Oregon that I want to display in my map sheets.
The imagery is surrounded by black pixels. I was able to remove the black pixels by checking the "Display Background Values: (R, G, B)" and setting it to 0, 0, 0.
But there is still (what I'm describing) as "noise" left around the pixels that I cannot get rid of (see image below).

I tried clipping the raster and that did not work.
I searched in several forums and couldn't find a resolution.
Any Ideas/suggestions?


texts: Generate polygon shaped as a circle if projected using Web Mercator, in Python?

texts:
I need to solve problem very similar to Generating polygon representing rough 100km circle around latitude/longitude point using Python? but subtly different.
I need a shape that will be a circle once projected in Web Mercator - not in reality.
Working answer ( https://gis.stackexchange.com/a/268251/45061 ) is not generalizable, as it relies on https://github.com/jwass/geog that will not allow doing what I am trying to do.


texts: I am interested in looping through a directory and listing only the red bands (Band 1) from 4-band NAIP imagery.  From there I would like to mosaic all of the band 1 layers into a new single band raster image using Mosaic to New Raster.  Is there a way to specify band 1 using arcpy.ListRasters() or should I find a new approach?


texts:
I deleted the layer heading thinking that I didn't need it. But after doing some refactoring, it turns out that I do need it but I can't figure out any way to add it back.

How do I add it back?

I'm using ArcMap 10.6.


texts: UpdateCursor stops after the first iteration

texts:
My arcpy script includes two updateCursor Objects. The first iteration of the "outer" coursor works fine, but the programm stops when it comes to the second iteration at updateRow(). Has somebody an idea, where the problem could be?

Here is the code snippet:

arcpy.env.workspace = "I:/GIS/VoGIS_Projekte/Wanderwege/Basisdaten/Basisdaten_Sommer/wandern_SK.gdb"

# Use row object to get and set field values
cursor = arcpy.UpdateCursor("wanderwege/wander_wege_3d")

# Iterate through rows and update values
row_i = 1
wa_j = 1
for row in cursor:
    print row_i
    row.WEG_LAENGE = '999'
    cursor.updateRow(row)               # The 2nd iteration stops here...
    wwege = arcpy.UpdateCursor("wanderwege/wander_wege")

    for wa in wwege:
        print wa_j
        wa.WEG_KATEGO = "Kaps999"
        wwege.updateRow(wa)                 # save
        wa_j = wa_j + 1
    del wwege
    del wa

    row_i = row_i + 1

del row
del cursor

I am working with:
ArcGIS 10.2
Python 2.7
Windows 7

P.S.: I know, that this code doesn't make much sense - I just wanted to reproduce the error of a quite complicated script.


texts:
I spend a comparison of functions of spatial analysis and Oraсle postgresql. There are two tables with the same coordinate system - tech_bridge.

oracle query

select sum(SDO_GEOM.SDO_AREA(t.gdo_geometry,1)) from geo.tech_bridge t

result : 167615.3831

postgresql query

select sum(ST_Area(t.gdo_geometry, true)) from geo.tech_bridge t

result : 167508.4431

why the results are different?


texts:
I want to run a python code that, if 4 fields are empty, insert in the first of them a text. The fields are all designed to contain text. My code is:

import arcpy
import xlrd
import numpy
import pdb
arcpy.env.workspace = r"Directory"
arcpy.env.overwriteOutput = True
shape = r"C:\Directory\Shapefiles.shp"
fields = ['A', 'B', 'C', 'D']
cursor = arcpy.da.UpdateCursor(shape,fields)
for row in cursor:
    if row[0] == None:
        if row[1]== None:
            if row[2]== None:
                if row[3]== None:
                    row[0] = 'Not_Officially_Recorded'
    cursor.updateRow(row)
arcpy.FeatureClassToShapefile_conversion(shape, r'Anotherdirectory')
print('Done')

Anyway, the code runs suspiciously quickly and after I look at it, I see that nothing has been done. So what may the error be?


texts:
I am making a WPS call in OpenLayers which works successfully in Chrome and Firefox but not in IE9. The response from IE tells me

  Process failed during execution\nCould not convert  to target type
  com.vividsolutions.jts.geom.Geometry

The same code returns intersecting polygons in both other browsers so I think it's down to how IE handles my GML/XML. Can anyone help out?

Code as follows:

function wpsExecute(featureA, featureB, identifier) {

    var polygonElement;
    if (BrowserInfo.browser == "Chrome") {
        polygonElement = "Polygon"
    } else {
        polygonElement = "gml:Polygon"
    }

    // set up the process specific parameters
    var options = {
        // set the wps process
        identifier : identifier,
        dataInputs : [ {
            identifier : "a",
            data : {
                complexData : {
                    attributes : {
                        mimeType : "text/xml",
                        subtype : "gml/3.1.1"
                    },
                    // get the gml:Polygon node from the xml DOM
                    value : featureA.getElementsByTagName(polygonElement)[0]
                }
            }
        }, {
            identifier : "b",
            data : {
                complexData : {
                    attributes : {
                        mimeType : "text/xml",
                        subtype : "gml/3.1.1"
                    },
                    // get the gml:Polygon node from the xml DOM
                    value : featureB.getElementsByTagName(polygonElement)[0]
                }
            }
        } ],
        responseForm : {}
    };

    var format = new OpenLayers.Format.WPSExecute();
    // console.log(format.write(options));

    try {

        OpenLayers.Request.POST({
            url : WPS_URL,
            data : format.write(options),
            async : false,
            success : parseExecute,
        });

WHERE featureA and featureB are as follows:

var gmlOptions = {
    featureType : "feature",
    featureNS : "http://example.com/feature"
};

// new gml 3 format
var gml3 = new OpenLayers.Format.GML.v3(gmlOptions);

compareFeatureGml = gml3.write(compareFeature);

var xmlFormat = new OpenLayers.Format.XML();
var featureB = xmlFormat.read(compareFeatureGml);


texts: Retrieve WMS service from GeoWebCache with OpenLayers

texts:
Ok guys, here we go !
For my needs, I basically added a standalone version of GeoWebCache.
After some configuration, the GeoWebCache retrieve the data from GeoServer using Capabilities.

Now, in my OpenLayers code, I would like to question GeoWebCache to retrieve my data, and not GeoServer (if it is possible).

For the moment, I have something like this (working) :

    var fond_carto = new OpenLayers.Layer.WMS(
        "Fond cartographique 2",
        "http://179.172.277.77:8080/geoserver/wms",
        {
            layers: 'departement',
            TILED: 'true'
        }

I tried stuff like that (not working) :

    var fond_carto = new OpenLayers.Layer.WMS(
        "Fond cartographique 2",
        "http://179.172.277.77:8080/geowebcache/sevice/wms",
        {
            layers: 'departement',
            TILED: 'true'
        }

Is there a way to accomplish that ?
After some reading, I assume I'm a bit lost with this standalone version of GeoWebCache. It is written it works as a proxy between client (OpenLayers) and server (GeoServer). But I read that the TILED: 'true' option means that GeoServer will send the layers for tiling to GeoWebCache. As I desactivated the integrated version of GeoWebCache in GeoServer, and as my WMS request are still tiled, I assume this is true. But then it's not really a proxy.

Could it works as follow:

openlayers => standalone geowebcache => geoserver => standalone geowebcache => openlayers

Or it has to work that way :

openlayers => geoserver => standalone geowebcache => geoserver => openlayers

Thanks for any tips or any clarifications :-) !

EDIT

After some more configuration in geowebcache.xml, I still have the error but I can add some details to help.
This is the URL of the request that is not displaying (pink tiles) :

http://179.172.277.77:8088/geowebcache/sevice/wms?LAYERS=departement&TILED=true&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image%2Fjpeg&SRS=EPSG%3A4326&BBOX=2.8125,48.515625,3.515625,49.21875&WIDTH=256&HEIGHT=256:8080/geowebcache/sevice/wms?LAYERS=departement&TILED=true&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image%2Fjpeg&SRS=EPSG%3A4326&BBOX=2.8125,48.515625,3.515625,49.21875&WIDTH=256&HEIGHT=256 :

and when I try to copy/paste this url, I got a GWC error :

400: Thread 35 Unknown layer departement. Check the logfiles, it may not have loaded properly.

I don't get why I got this error. Here's my geowebcache.xml description of the layer :

<wmsLayer>
  <name>departement</name>
  <mimeFormats>
    <string>image/gif</string>
    <string>image/jpeg</string>
    <string>image/png</string>
    <string>image/png8</string>
  </mimeFormats>
  <gridSubsets>
    <gridSubset>
      <gridSetName>EPSG:4326</gridSetName>
    </gridSubset>
  </gridSubsets>
  <wmsUrl>
    <string>http://179.172.277.77:8080/geoserver/wms</string>
  </wmsUrl>
  <wmsLayers>aaz:departement</wmsLayers>
</wmsLayer>


texts:
I need to implement a WFS Server using Geo tools, I have an application developed in Java using a spatial database.  I wish to develop  a WFS service using GeoTools. As a first step, I am implementing the "GetCapabilities" functionality  and generating information in GML,  both GET and POST queries.

My idea is to use the Geotools classes, which will allow me to implement the following parts: Identification Service, Service Provider, Operation Metadata, etc., and then to generate a GML document to send to the client.

In using Geotools, the Geotools class started with GetCapabilitiesResponse (org.geotools.data.ows) and Class WFSCapabilities (org.geotools.data.wfs.v1_0_0).
Could someone please help me with a reference, tutorial or a sample code, to know how to use it?

Patrick Morales


texts: Finding geographic coordinates of shapefile using ArcPy with ArcMap

texts:
I want my Python script to check if the first coordinate in a GeoJSON polygon is within a polygon in ArcMap, and if so draw the polygon to a featureclass. My issue is that the GeoJSON file is using GCS while the coordinates I get from my shapefile are in PCS.

desc = arcpy.Describe(boundary)
xmin = desc.extent.XMin
ymin = desc.extent.YMin

if xmin < Geopoint and ymin < Geopoint:
    with arcpy.da.InsertCursor(feature_class, ['SHAPE@']) as cursor:
        cursor.insertRow([lst])
else:
    pass

How can I get the x and y min or max of my boundary shape in PCS so that I can compare it to my GeoJSON point which is in GCS?

I feel like there must be a ArcPy function but I cannot seem to find it.

I am also aware that there may be simple ways to do this in ArcGIS Pro, however my organization does not use ArcGIS Pro.


texts: What is the fastest way to iterate features for geometry operations?

texts:
There was an old GIS soft "ArcView" ; scripts in avenu language which worked much faster than my Add-ins now (made in Visual Studio in С#). So I would like understand what did I do wrong? Why is the modern technology less effective than the old one?

The code calculates intersection area of every combination of polygons in two featureclasses.
Featureclass fc1 contains 100 features, FC fc2 contains 1000 features, therefore there are 100,000 combinations.

My first try was "FOR" cycle (without using cursors) this way:

IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
IActiveView activeView = mxdoc.ActiveView;
IMap map = activeView as IMap;

ILayer layer1 = map.Layer[0];
ILayer layer2 = map.Layer[1];

IFeatureClass fclayer1 = (layer1 as IFeatureLayer2).FeatureClass;
IFeatureClass fclayer2 = (layer2 as IFeatureLayer2).FeatureClass;

IFeature f1;
IFeature f2;
IFeature fINTERS;

IGeometry g1;
IGeometry g2;
IGeometry gINTERS;

IArea arINTERS;
double areaINTERS;

ITopologicalOperator topoINTERS;

for (int i = 0; i < fclayer1.FeatureCount(null); i++)
{
    f1 = fclayer1.GetFeature(i);
    g1 = f1.Shape as IGeometry;
    topoINTERS = g1 as ITopologicalOperator;

    for (int j = 0; j < fclayer2.FeatureCount(null); j++)
    {
        f2 = fclayer2.GetFeature(j);
        g2 = f2.Shape as IGeometry;
        gINTERS = topoINTERS.Intersect(g2, esriGeometryDimension.esriGeometry2Dimension);
        fINTERS = fclayer1.CreateFeature();
        fINTERS.Shape = gINTERS;
        arINTERS = fINTERS.Shape as IArea;
        areaINTERS = arINTERS.Area;
        fINTERS.Delete();
    }
}

It takes about 90 minutes to do this calculation!

My second try was "WHILE" cycle (using cursors) this way:

IFeatureCursor fcur1 = fclayer1.Search(null,true);
IFeatureCursor fcur2;

IFeature f1 = fcur1.NextFeature();
IFeature f2;
IFeature fPEREC;

IGeometry g1;
IGeometry g2;
IGeometry gINTERS;

IArea arINTERS;
double areaINTERS;

ITopologicalOperator topoINTERS;

while (f1 != null)
{
    g1 = f1.Shape as IGeometry;
    topoINTERS = g1 as ITopologicalOperator;

    //update cursor
    fcur2 = fclayer2.Search(null, true);
    f2 = fcur2.NextFeature();

    while (f2 != null)
    {
        g2 = f2.Shape as IGeometry;
        gINTERS = topoINTERS.Intersect(g2, esriGeometryDimension.esriGeometry2Dimension);
        fPEREC = fclayer1.CreateFeature();
        fPEREC.Shape = gINTERS;
        arINTERS = fPEREC.Shape as IArea;
        areaINTERS = arINTERS.Area;
        fPEREC.Delete();

        f2 = fcur2.NextFeature();
    }
    f1 = fcur1.NextFeature();
}

This code spends about 2 minutes to do this, it's too long for my purpose too.
Script in avenue (ArcView) spends about 1 second on this task (script uses "FOREACH" cycle, but in C# is not allowed to iterate the features within featureclass this way).

Is there any idea how to optimize this code to make it work faster?

Using: C#, ArcGIS 10.2.2


texts: Conditional statement in R based on multiple rasters

texts:
I have:

raster_1: an integer raster with non-sequential values i (101, 130, 240, 323)
i: a watershed discharging X of sediment in the environment
raster_2: represents the water quality for each watershed i with values ranging from 0 to 100,000 - named plume_i

I have ~400 rasters representing the water quality of 400 watersheds.
I would like to build a loop to set the values in raster_1 to NA for each zone (defined by the values: 101, 130, 240, 323 in raster_1) when raster_2 for the corresponding value (e.g., plume_101) is 0, using R.
For instance,

When raster_1 = 101 & plume_101 = 0, assign NA to raster_1 otherwise return raster_1 values
When raster_1 = 130 & plume 130 = 0, assign NA to raster_1 otherwise return raster_1 values
When raster_1 = 240 & plume 240 = 0, assign NA to raster_1 otherwise return raster_1 values
and so forth...

How can I do this?


texts: QGIS (2.0.1): Field calculator for length not working - values way too high

texts:
Same problem as another thread: Field Calculator returning incorrect length and area values

Trying to calculate length of a line shapefile, and field calculator is returning values that are way too high (by factor of 1000 or more, like the thread above). Both the shapefile and the project are projected in NAD83 / North Carolina (ftUS) - so they are consistent. The preview shows lengths that appear correct/reasonable, but then the calculated value that is returned is way off.

See attached image - shows that output preview looks right, but then the calculated values (to the right) are way higher and don't match the output preview.

Any suggestions?


texts: I would like to select the largest number in a column of an attribute table by "select by attributes" in QGIS to select the polygon with the largest area. How would be the expression for that? I do not want to sort the column and select the polygon per hand, because I want to use it in the modeler.


texts: How to implement Shortest Path in Openlayers 3 and Geoserver 2.8.2?

texts:
I've managed to get all the SQL queries necessary to get the shortest path with Dijkstra using pgRouting, but now I want to implement them in my website using OpenLayers 3 and GeoServer 2.8.2

The first thing that got in my mind was to use the SQL Views but I don't know whether it's the best solution or if it is possible.

Then I thought about using WPS but I'm not used to it.

So if anyone have already faced the same problem could you help please?


texts: QGIS layout manager has several A4 layouts: export them all to "the same pdf"

texts:
My QGIS layout manager has several A4 layouts.
I would like to export them all to "the same pdf".

I think that by going through the python console it should be possible.
for now I can only print one page....
Here is my code:

manager = QgsProject.instance().layoutManager()
exporter = QgsLayoutExporter(layout)
for layout in manager.printLayouts():
      print(layout.name())
      exporter.exportToPdf("F:/44.pdf",QgsLayoutExporter.PdfExportSettings())

Here is the return: 0, 1, 0, 2

Only page 2 is exported to PDF


texts: Setting QgsExpression for QgsPalLayerSettings in QGIS

texts:
How do you apply a QgsExpression to a Label ?

I can do it in QGIS with this simple expression
concat('N:' ,round(Value,0),'\n','Yield:' ,round(Value 2,0))

N: 65

Yield: 40

How do I do it in PyQGIS ? I see a getLabelExpression but no setLabelExpression

This works but I want but only for the field name

 crap = QgsVectorLayer("D:/data.shp", "data", "ogr")
 QgsMapLayerRegistry.instance().addMapLayer(crap)
 palyr = QgsPalLayerSettings()
 palyr.readFromLayer(crap)
 palyr.enabled = True
 palyr.fieldName = 'Name'
 palyr.placement= QgsPalLayerSettings.OverPoint
 palyr.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'8','')
 palyr.writeToLayer(crap)


texts: QGIS Geometry $area Field calculator function compared with Cadcorp SIS

texts:
I want to know what QGIS uses as a measurement when using the Geometry and $area function. As in is, the area calculated Sq meters?

I am getting two different answers when using Cadcorp SIS and QGIS:

Cadcorp SIS: 132068457.409951

QGIS: 131747637.35

My session is in British National Grid EPSG:27700.


texts: turn osm2po vertex table restrictions column into turn cost for shooting star pgrouting?

texts:
i found this thread

Does OSM2PO take into consideration turn restrictions?

but im kind of new to sql and database so me, making my own query or sp to achieve the above suggestion by myself is kind of far from reality. can somebody help me please? thanks!


texts: How to extend an OSM or raster area in QGIS (into 'No Data' area)?

texts:
I need to do the opposite of clipping down a raster (or SHP file) in QGIS (or Gdal, Python etc - even R if someone can provide the script).

I have a dataset from Open StreetMap (*.OSM) which includes the south west corner of a country. I have also converted this to a SHP file or Raster format).

To match with my other Maxent datasets I need to expand the OSM data south and west so I am only increasing the 'no data' areas at the left and bottom of the files.

I have looked at previous questions but not found a clear answer for this. Is there a way to do this in QGIS?


texts:
I have a geoserver installed on my local machine and on it i have published the map which is connected to Postegresql.

Now i can view a map on geoserver using the openlayers viewer (url-"http://localhost:8081/geoserver/test/wms?service=WMS&version=1.1.0&request=GetMap&layers=test:city_limit&styles=&bbox=-114.315796136219,50.842822033593,-113.859905274451,51.2124253170889&width=768&height=622&srs=EPSG:4326&format=application/openlayers")"

Now i have code api of leaflet to customise a map which is hosted on geoserver.

Here is a HTML code:

    <!DOCTYPE html>

    <title> Leaflet Demo </title>

   <link rel="stylesheet" href="https://npmcdn.com/leaflet@1.0.0-rc.3/dist/leaflet.css" />

    <form id="form1" >

    var map = L.map('map').setView([55.67, 12.60], 11);

    var mywms = L.tileLayer.wms("http://localhost:8081/geoserver/test/wms", {
    layers: 'k101:theme-startkort',
    format: 'image/png',
    transparent: true,
    version: '1.1.0',
    attribution: "myattribution"
    });
    mywms.addTo(map);

    </form>

But when compiling on visual studio ,browser displays :
a code itself

Any idea, what to do?


texts: Does anyone already parallelized mosaic function for rasters in R?

texts:
Does anyone already paralleled mosaic function for rasters in R?

It seems that snow still doesn´t deal with this on an automatic way.

The idea is something like this:

r1<-stack("fit__test23KLR.img")
r2<-stack("fit__test23KKR.img")
r3<-stack("fit__test23KKS.img")
r4<-stack("fit__test23KLS.img")

ff <- function(x1, x2, x3, x4) mosaic(x1, x2, x3, x4, fun="mean")

beginCluster(35)
cl <-getCluster()
m1 <- clusterR(r1, r2, r3, r4, fun = ff, export="mosaic")
endCluster()

Each stack r1, r2.... are:

dimensions  : 2516, 412, 1036592, 93  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
range between: -2.5 to 2.5


texts: Why are my modifications to UniqueValuesSymbology classLabels not following through to PDF export?

texts:
I have a vector layer that I add to a data frame. I modify the unique values symbology class labels. I then export the map document as a PDF and the labels in the PDF do not reflect the modifications.

Here is the code and the resulting print statements:

#Dissolve zones and get a list of acres
arcpy.Dissolve_management(pathScratchGDB + "fcZones_Proj", pathScratchGDB + "fcSurface_Temp",
    dissolve_field="ph1", statistics_fields="acres SUM", multi_part="MULTI_PART", unsplit_lines="DISSOLVE_LINES")
lsAcres_Temp = list(row.getValue("SUM_acres") for row in arcpy.SearchCursor(pathScratchGDB + "fcSurface_Temp", fields="SUM_acres"))

#Format symbology
df = arcpy.mapping.ListDataFrames(mxdSurface, "df1")[0]
df.spatialReference = srCover
legend = arcpy.mapping.ListLayoutElements(mxdSurface, "LEGEND_ELEMENT", "df1_legend")[0]
legend.autoAdd = True
lyrTemp = arcpy.mapping.Layer(pathScratchGDB + "fcZones_Proj")
arcpy.ApplySymbologyFromLayer_management(lyrTemp, pathLayerTemplates + "/zone_surface.lyr")
lyrTemp.symbology.valueField = "ph1"
lyrTemp.symbology.addAllValues()

#Label zones
lblClass = lyrTemp.labelClasses[0]
lblClass.expression = "\"<FNT size = '8'>\" + " + "[zone_num]" + " + \"</FNT>\""
lyrTemp.showLabels = True
lblClass.showClassLabels = True

#Add Layer
arcpy.mapping.AddLayer(df, lyrTemp, "TOP")
#Update legend style
for lyr in legend.listLegendItemLayers():
    legend.updateItem(lyr, styleItem)

#Modify legend to include acres
print(lyrTemp.symbology.classLabels) #[u'7.1', u'7.2', u'7.4']
temp = lyrTemp.symbology.classLabels
iNumLabels = len(temp)
for i in range(0,iNumLabels):
    temp3 = str(temp[i]) + " (" + str(round(lsAcres_Temp[i],1)) + " ac)"
    temp[i] = str(temp3)
lyrTemp.symbology.classLabels = temp
print(temp) #['7.1 (33.0 ac)', '7.2 (43.4 ac)', '7.4 (69.4 ac)']
print(lyrTemp.symbology.classLabels) #[u'7.1 (33.0 ac)', u'7.2 (43.4 ac)', u'7.4 (69.4 ac)']
arcpy.RefreshActiveView()

#Export PDF
arcpy.mapping.ExportToPDF(mxdSurface, "../Scratch/temp.pdf")

And here is an image snipped from the PDF:

Why don't the modified class labels follow through to the PDF?

It might be worth noting that I've successfully exported a PDF with the modifications when I do everything in the ArcMap python window.


texts:
I'm working on an old table where polygons are stored as points.
Now, I want to migrate these points to PostGIS and store them as polygons (or multi-polygons).
The table structure is as follows:
loc_id; loc_shape_id; geom;

Here the points for loc_id 29930:

If I now apply a simple ST_MakeLine like the following:

select
st_makeline(geom) as line, loc_id
from (select * from test_regionen_11
order by loc_id, loc_shape_id )
as ordered_points
where loc_id = 29930
group by loc_id

I get this messed up linestring:

Well, does anyone has an idea how I could create multi-linestrings that represent the island in this case?

There is no column or value to distinguish them.


texts: ArcGIS 10.1: "ERROR 000223: Cannot extract event table properties" when iterating through multiple files/folder

texts:
I am trying to convert ~270k csv-files with lat/long-data to layerfiles. The files are split into folders by year, i.e. ~30 folders with ~9k csv-files in each.

I have written a script to iterate through the files and it works well for test-samples of 100+ split into a couple of different folders. However, when I try to run it on the full set of files, it will either:

work well for the first folder (~9k files) and then when moving on to the second folder it will give the following error for all subsequent files:
arcgisscripting.ExecuteError: ERROR 000223: Cannot extract event table properties
Failed to execute (MakeXYEventLayer).
on my other computer it will usually stop after an arbitrary amount of file-conversions before it finishes the first folder (typically at +-8.5k of the conversions) and then give the same error for all subsequent files

The error seems to be similar to this: http://forums.arcgis.com/threads/80161-Make-XY-Event-Layer-Error ,but my error pops up after an arbitrary amount of files has been processed.

I've included my code for reference below. Note that I have also tried to split the list of files per folder into chunks of 100 and defined the layercreation as a function that is being called for each chunk, but the error still persisted in the same way. I suspect that it might be some setting in ArcGIS or on my computers that causes it, but I don't know what. Does anybody have any idea where I might start to solve this?

#setting paths for the project
path="C:/project"
savepath = path + "/outputlayers"
csvpath = path+ "/csv"

#lists the directories in the csv directory
yearlist = os.listdir(csvpath)

#checks the directory of each year:
for year in yearlist:
    indir = csvpath + "/" + year
    hourlist = os.listdir(indir)

        #runs through each file and converts it to a .lyr file in the outputdirectory
        for infile in hourlist:
            try:
                # Set the local variables
                in_Table = indir + "/" + infile
                x_coords = "longitudex"
                y_coords = "latitudey"
                z_coords = ""
                out_Layer = infile[:-4]
                saved_Layer = savepath + "/" + year + "/" + infile[:-4] +'.lyr'

                # Set the spatial reference
                spRef = arcpy.SpatialReference(4301)

                # Make the XY event layer...
                arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer, spRef, z_coords)

                # Print the total rows
                print arcpy.GetCount_management(out_Layer)

                # Save to a layer file
                arcpy.SaveToLayerFile_management(out_Layer, saved_Layer)

            except:
                # If an error occurred print the message to the screen
                print arcpy.GetMessages()


texts: Creating separate line feature from multi-line feature in QGIS?

texts:
I am working with road network data and I have a large amount of lines that are not connected but are classified as the same feature.  I want to select one or more of those lines and create a new feature in order to rename and connect it with the correct road feature.

Is there an easy way to do this?

I've found a workaround that requires saving the selected feature as a layer, running 'single parts to multiparts', selecting the features I want, creating another layer, then merging that layer with the original (which creates a 4th layer).  Not exactly efficient.


texts: Add minutes contained in a column to a date

texts:
In CartoDB I want to add minutes contained in a column of data type number to a date.
e.g.

SELECT my_initial_date + 'xx minute'::INTERVAL as my_date FROM table

What is working is:

SELECT my_initial_date + '7 minute'::INTERVAL as my_date FROM table

This adds 7 minutes to each my_initial_date.
What I want to do is the same but with a variable number of minutes stored in another column, e.g.:

SELECT my_initial_date + 'column_name minute'::INTERVAL as my_date FROM table

This doesn't work in CartoDB! How to do this?


texts: Convert *.asc OSGB36 files to WGS84 GeoTIFF using gdal/ogr2ogr

texts:
I am trying to obtain elevation readings using gdallocationinfo (as per other datasets I have) from a set of .asc files containing elevation data (plus an associated .prj file), using WGS84, but I can't seem to find a way to do this in a way that uses the correct projection.

gdal_merge.py -o ../osTerrain50.tiff *.asc

gives me a GeoTIFF.

Opening this in QGIS shows that the data appears correctly, but it does ask me for the projection, so I specify EPSG:27700 (which is OSGB36).

Indeed gdalinfo does show that the .tiff file doesn't end up containing the projection:

> gdalinfo osTerrain50.tiff

Driver: GTiff/GeoTIFF
Files: osTerrain50.tiff
Size is 13200, 24600
Coordinate System is `'
Origin = (0.000000000000000,1230000.000000000000000)
Pixel Size = (50.000000000000000,-50.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (       0.000, 1230000.000)
Lower Left  (   0.0000000,   0.0000000)
Upper Right (  660000.000, 1230000.000)
Lower Right (  660000.000,       0.000)
Center      (  330000.000,  615000.000)
Band 1 Block=13200x1 Type=Float32, ColorInterp=Gray

What do I need to do reproject this to WGS84, into some format (TIFF/shp/whatever) that gdallocationinfo can read?

NB I need to do this via the command-line, not via GUI programs.


texts: Zooming to extent of filtered vector layer using PyQGIS?

texts:
I have a vector layer that I need to filter and zoom to extent of the layer whilst the filter is applied. Most of the examples I see only zoom to selected features. I have managed to filter using

p_id = 31
filter_rule = ' id = %s ' % p_id
qExp = QgsExpression(filter_rule)
qReq = QgsFeatureRequest(qExp)
features = layer.getFeatures(qReq)

How can I zoom to the extent to only show the feature which corresponds to the filter?


texts: Batch-enable archiving on multiple datasets

texts:
The ArcGIS Desktop 10.7
 Enabling Archiving page says:

  To enable archiving on multiple datasets at once, you can use the
  Enable archiving geoprocessing tool.

When I open the Enable Archiving tool in ArcGIS Desktop 10.7.1, I don't see a way to select multiple datasets:

How can I enable archiving on multiple datasets at once using out-of-the-box functionality?


texts: How to compile in Fortran in order to use an add-on in GRASS?

texts:
I need to use the r.inund.fluv add-on in GRASS. I installed the gfortran compiler in order to use the script. The problem is, I do not know how to compile the files in gfortran. Im using windows

Here are the files for r.inund.fluv.

I emailed the author of the script, Roberto Marzocchi and he instructed me to do the following and I was stuck at step 4.:

Try the following way:

1) save the file r.inund.fluv within the folder named scripts within the grass main folder

2) create a file r.inund.fluv.bat in the folder named bin and copy the following row:

@"%GRASS_SH%" -c '"%GISBASE%/scripts/r.inund.fluv" %*'

3) save the file r.inund.fluv.html within the folder named doc

4) create a folder named fortran_code within your etc folder and save the four fortran code 2d_path.f90, etc etc. then try to compile the 4 fortran codes using for the binary file the same name of the source code without any extension  e.g (2d_path.f90 become 2d_path)


texts: Zooming to particular feature among selected features by feature id in PyQGIS?

texts:
Using QGIS 2.6 and PyQGIS, I want to zoom to features one by one on click among selected features list and change the colour on zoom without unselecting the other features.

canvas=iface.mapCanvas()
canvas.selectToFeaturesIds(layer,id)   # not supported in QGIS 2.6

or

canvas.panToFeatureIds(layer,id)     # not supported in QGIS 2.6

I have enlisted the selected features in Qtablewidget and want to zoom to them one by one by mouse click without unselecting others


texts: Changing the target field width using Fieldmappings. ArcGis, Merge,

texts:
I need to merge several shapefiles. All of them has the same attribute/field structure, but the width is not the same for the attributes with the same name.
I tried to use Fieldmappings to set a wide enough destination field. But is not working. It doesn't create a long (50 characters) enough, so the error persistes.
Error:  Error 001156: Failed on input 26788, could not write value 'looong value here'. I know that the long strings are the problem. I've erased that polygon ID, and the error is gone.
However, I need to have this process automated, solving it manually is not the solution.
What is wrong with my code/Fieldmappings?

folder_combined_rasters        = r"C:\Projects\Soils\FECv3Out.gdb"
MapSeries   = []

# loading all feature classes into shapefiles memory list
arcpy.env.workspace = folder_combined_rasters
shapefiles = arcpy.ListFeatureClasses()
shapefiles.sort()

# Adding only the shapefiles that ends in 'Series' into MapSeries
for j,county in enumerate(shapefiles):
    if county[3:] == 'Series':
        MapSeries.append(county)

# defining fieldmaps
fldmps      = arcpy.FieldMappings()
fm_series   = arcpy.FieldMap()

# Adding only the needed fields to the fieldmap output
for field in arcpy.ListFields("SA_Series"):
    field.lenght = 60
    fm_series.addInputField(fm_series, field.name)

fldmps.addFieldMap(fm_series)

seriename = "Mine_Series"

arcpy.Merge_management(MapSeries,   seriename , fldmps)


texts: Using Field Mapping with TableToTable_conversion in arcpy

texts:
I am importing csv files into a gdb, and need to change a couple field types during the import. I am trying to use the field mapping in TableToTable_conversion to change the field type as the table gets imported.

I am confused about how to identify the field that needs to be updated. When setting up the field mapping I am using the fm.addInputField line, which requires a field name. But because the gdb table does not exist yet as it has not yet been imported, I can't define the field to input.

Should the approach be to use the csv rows to define which field I want to change, or am I missing something completely? Is there an easier way to change a field type upon import rather than using Field Mappings?


texts: Modify 'Codepoint with Polygons' to add new fields

texts:
I'm using OS Codepoint With Polygons data set which provides shapefiles of each Postcode area within the UK. Within each postcode area are boundaries of individual postcodes.

I'm trying to merge these layers together to make another shapefile of our business boundaries that are made up from groups of postcodes.

From reading on here, it seems that I need to add a field to the Dbf file to name the business area that I want to 'dissolve' to and fill that field with the business area attributes.

Am I barking up the right tree here?


texts: Buffering based on topography in ArcGIS Desktop?

texts:
I need to do a buffer based on the real terrain.
Like 500 meters from a polyline on all sides.
I have DEM (SRTM) of the terrain.

Is there a way to do that?


texts: Remove gaps between two polygons using ArcGIS for Desktop?

texts:
I have already found some ideas how to solve this problem but none of them seem to fit perfectly.

I need to remove the gaps between a shapefile and a Feature Class in our Database so that I can import the shapefile without topology differences. Here is an example:

The shapefile polygons are mostly parcels of the bigger one. The problem mainly occurs when there are curves. So I need to adjust that curved border without changing the upper border. And this for like 2000 polygons.

I found out about 3 possibilities:

Snap tool -> wont change anything
integrate tool -> looked pretty fine in the beginning but it also changes the inner borders
trace tool -> tracing the borders manually is a way but before I start this with 2000 polygons I would like to ask my question here.


texts: How to calculate area correctly in QGIS?

texts: I try to measure area in QGIS, but encounter some problems. The results seem to be quite inexact. When testing it with some administrative units (polygons from OpenStreetMap) and compare the results with the known real surface, results in some cases are quite close, while in other cases differ up to 20% or more from the "real" surface. So there seems to be a problem with my measurement. I saw some answers in this forum connected to projection choice. I couldn't understand  which projection (the name I have to choose) is adequate for area measurements?


texts: Generalizing Polygons to MultiPolygons in GeoDjango?

texts:
I set up a model with models.PolygonField in geodjango, using postgres as database. I try to import shp into postgres. The problem is, shp (compiled with QGIS) mixes polygon and multipolygon, hence it always fails to do the export because of the constraint check enforce_geotype.

Is there a way to clear the constraint, so as to store both multipolygon and polygon type data?


texts:
I'm very new to the PostGIS and QGIS. I'm a student and I am doing thesis work on map matching.
I have segment shapefile of two network one OSM network and other one TMC network. And I have point shapefile of two networks.
Methodology to be done;

I need to find the closest point of OSM points to the TMC point. I have find using closet point buy using distance matrix. The distance matrix attribute table contains  inputid as TMC points id and target id as OSM points id and also I also have distance between them.

In OSM network I don't have start and end point of the segments. How can I get start and end points of the segments?

And the my goal is to map TMC SINGLE LINK to OSM PATH.

For each single link between TMC POINTS, I have to find (MAPPING) which is the correspondent path (SEQUENCE OF LINKS) on OSM NETWORK, between the correspondent points OF TMC NETWORK on the OSM NETWORK. On OSM links should be attached the attribute of the link ID OF TMC  ELEMENT.
For example, suppose to have the following situation:
TMC NETWORK:     POINT 1: ID =500; POINT 2: ID=501; LINK between POINT 1 and POINT 2: ID = 7500
We found that on OSM NETWORK, the correspondent POINT of POINT 1 with ID =500 on TMC, is the OSM POINT 1 with ID = 12349876
We found that on OSM NETWORK, the correspondent POINT of POINT 2 with ID =501 on TMC, is the OSM POINT 2 with ID = 55532481
Now, we have to calculate, on OSM NETWORK,  which is the shortest path between the POINT 1 (ID = 12349876 ) and the POINT 2 (ID = 55532481 ) on OSM NETWORK.
Suppose that we found the  SEQUENCE OF LINKS (78934, 78564, 12789, 35785).
At the end, we can assume that the arch with ID = 7500 on TMC NETWORK, corresponds to the SEQUENCE of ARCHES with ID = ( 78934, 78564, 12789, 35785). So we can add to OSM ARCHES ( ID = 78934, ID = 78564, ID = 12789, ID = 35785) the  attribute "TMC ARC ID" with the value: 7500.

Can you tell me what are procedure to be done for points 2 and 3  by using QGIS and PostGIS, pg routing.


texts: Calculating flowpath length to stream within catchment using ArcGIS Desktop?

texts:
ArcGIS 10.3 provides a hydrology tool "flow length" which calculates the length of flow paths within a basin to the outlet of the basin. I am trying to calculate flow path length and flow path gradient from a given pixel of interest to where it enters the stream network. Using TauDEM (which I do not have access to), they define this as the 'distance from a pixel of interest to the stream channel using the D8 flowpath direction'.

I have a 30m DEM, flow direction raster and river network. From these flow path values I need to calculate mean flowpath length and mean flowpath gradient for a catchment for use in a catchment classification study.

Using ArcGIS 10.3 for Desktop.


texts: Create Feature Class from another Feature Class using ArcObjects for .net application in C#?

texts:
I have a Feature Class and I want to create another Feature Class that is basically a subset of the first Feature Class. For example, I have a Feature Class containing public parks which I want to use to create a new Feature Class containing only certain parks.

Should I be using the IQueryFilter interface or is there a better/easier way to do this?

I am creating a .net application in C#.


texts: Calculate field python script returning NULL values

texts:
I have a script that I am using to create and calculate a new field based on the features spatial location. Basically, I have a geodatabase that contains parcels by county in florida, and I am trying to calculate which wind_zone they are in based on their location within a shapefile that contains wind zone polygons. The script runs fine, but the wind zone description values are not transferring over, instead the field is populated with nothing but NULL values.

The code I have so far:

import arcpy

# Set overwrite option
arcpy.env.overwriteOutput = True

# Define your workspace. This should be your database with parcels feature     class
arcpy.env.workspace = r"C:\Data\HAAS\parc_join\parcels_1.gdb"

# This can be done once for all iterations
arcpy.MakeFeatureLayer_management(r"C:\Data\HAAS\wind_zones\WZ_merge.shp",     "lyr_windZones")

for fc in arcpy.ListFeatureClasses():
    # KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)
    # Create FeatureLayers
    arcpy.MakeFeatureLayer_management(fc, "lyr_parcels")
    # Add a Windzone description field
    arcpy.AddField_management("lyr_parcels", "WZ_Desc1", "SHORT")

    # Create a search cursor for the WZs
    rows = arcpy.SearchCursor("lyr_windZones")
        for row in rows:
        # What you'll do is select each WZ one at a time, and then select all the parcels in that WZ and calculate the fields
        # NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
        arcpy.SelectLayerByAttribute_management("lyr_windZones", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
        arcpy.SelectLayerByLocation_management("lyr_Parcels", "HAVE_THEIR_CENTER_IN", "lyr_windZones", "", "NEW_SELECTION")
        arcpy.CalculateField_management("lyr_parcels", "WZ_Desc1", "'{0}'".format(str(row.getValue("DESCRIPT"))), "PYTHON_9.3", "")
        print "Finished processing " + fc +str(row.getValue("DESCRIPT"))

print "Script Complete"


texts: canvas vs webgl client library to draw map vector tiles

texts:
For a map application using pbf vector tiles, should I go with Tangrams which uses webgl to draw the tiles or some other library (e.g. Hoverboard)?

Tangrams looks very promising since it abstracts out the tile styling to a .yaml file, but it won't work on clients without webgl support, and the ones on which it will, I'm afraid it'll consume more resources than a library using canvas.

So on a related note, which would be better - using canvas or webgl to render a map?


texts: gdal python set projection of a raster not working

texts:
Gdal 1.9.2, ECW driver: ECW/ERDAS Compressed Wavelets (SDK 3.x)

I am having problems with the script that supposed to set projection of the ECW to EPSG:27700, i.e. the British_National_Grid.

Please could you have a look at the script and suggest what could be wrong?
I created two versions of teh script and neither of them work :(

I have checked and it does not seem to be the ECW file format/ECW sdk issue as the script does not updte the projection for Geotiffs either.

script 1:

import sys, os, os.path, string, array, math, shutil, glob, re, locale, math, traceback, time, datetime, csv
from osgeo import gdal, gdalconst, osr
from osgeo.gdalconst import *
import subprocess

driver = gdal.GetDriverByName("ECW")
driver = gdal.GetDriverByName("GTiff")

# set OSGB96
sr = osr.SpatialReference()
sr.SetProjection ("EPSG:27700")
sr_wkt = sr.ExportToWkt()
print sr_wkt

inDir = str(sys.argv[1])
for file in os.listdir(inDir):
    ds = gdal.Open(inDir + "\\" + file, gdal.GA_Update)
    if ds:
        print('Updating projection for ' + file)
        res = ds.SetProjection(sr_wkt)
        if res != 0:
            print('Setting projection failed ' + str(res))
        ds = None # save, close
    else:
        print('Could not open with GDAL: ' + file)

and script 2:

import sys, os, os.path, string, array, math, shutil, glob, re, locale, math, traceback, time, datetime, csv
from osgeo import gdal, gdalconst, osr
from osgeo.gdalconst import *
import subprocess

driver = gdal.GetDriverByName("ECW")
driver = gdal.GetDriverByName("GTiff")

ds_reference = gdal.Open('D:\\reference_geotif\\reference.tif', gdal.GA_ReadOnly)
prj_wkt = ds_reference.GetProjectionRef()
print prj_wkt

inDir = str(sys.argv[1])
for file in os.listdir(inDir):
    ds = gdal.Open(inDir + "\\" + file, gdal.GA_Update)
    if ds:
        print('Updating projection for ' + file)
        res = ds.SetProjection(prj_wkt)
        if res != 0:
            print('Setting projection failed ' + str(res))
        ds = None # save, close
    else:
        print('Could not open with GDAL: ' + file)


texts: SELECT assign_vertex_id('ways', 0.00001, 'the_geom', 'gid'); not work

texts:
I'm using using OSgeo 6.5 dvd and studing pgrouting workshop (link) for get shortest path. I created network topology call "pgrouting - workshop". There are three commands to execute in section 5.2 they are :

-- Add "source" and "target" column
ALTER TABLE ways ADD COLUMN "source" integer;
ALTER TABLE ways ADD COLUMN "target" integer;

-- Run topology function
SELECT assign_vertex_id('ways', 0.00001, 'the_geom', 'gid');

Add source and target was executed successfully. topology function give errors.
Heres the error :

NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 15 at EXECUTE statement
ERROR:  function startpoint(geometry) does not exist
LINE 1: SELECT gid AS id, StartPoint(the_geom) AS source, EndPoint(t...
                      ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
QUERY:  SELECT gid AS id, StartPoint(the_geom) AS source, EndPoint(the_geom) as target FROM ways
CONTEXT:  PL/pgSQL function "assign_vertex_id" line 27 at FOR over EXECUTE statement

How do i fix this error?

I loaded the legacy.sql and noted that error has little bit changed :

NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 15 at EXECUTE statement
ERROR:  function point_to_id(geometry, double precision) does not exist
LINE 1: SELECT point_to_id(setsrid(_r.source, srid), tolerance)
           ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
QUERY:  SELECT point_to_id(setsrid(_r.source, srid), tolerance)
CONTEXT:  PL/pgSQL function "assign_vertex_id" line 33 at assignment


texts: ArcPy not executing CopyFeatures_management to write feature class to ArcSDE geodatabase

texts:
I have some problems with my Python script. The script gets an JSON with different information about images. The information includes general information about images and also the gps position of them. Wih this information a featureclass should be generated that shows the location of the pictures as point geometry. The other information should be stored in the attribute table.

The script runs very well and is already deployed on ArcGIS Server 10.1. Now the problem is, that at the first start the featureclass has to be created. As the featureclass consists of about 92 attribute fields the process of creating the featureclass is done with the option "in memory" as this option guarentees a fast creation of the featureclass - it is at least faster then directly running arcpy.addfield_management() on a feature class that was created the other way.

The disadvantage is that I want to store the featureclass permanently in a ArcSDE geodatabase. For this I have to run one of the tools arcpy.featureclasstofeatureclass_conversion or arcpy.copyfeatures_management. If the script is trying to run on of these lines I get the error message 000210 that is telling me that the output cannot be generated.

The connection to the ArcSDE works fine, this cannot be the problem. I have full write and read access. Does anybody know what the error message is about and what could be one solution?

The script works fine if I am not using the "in memory" option for the feature class. So this error must have his origin in the "in memory" option. The log file of the Database gives only an general error message:

  Oracle : ORA-3137[3149]


texts: QGIS Adding categorial color to points from CSV

texts:
I'm new to QGIS and couldn't find any answers to this question I had.

I'm importing data from a CSV file into QGIS using Layer > Add Layer > Add Delimited Text Layer. Here are a few lines of sample data:

  1,0,-72.27141,42.564782,ATHOL,15 1,1,-72.27141,42.573775,ATHOL,15
  1,7,-72.25924,42.573775,ATHOL,15 1,8,-72.25924,42.582768,ATHOL,15
  1,12,-72.24707,42.564782,ATHOL,15 1,14,-72.24707,42.582768,ATHOL,15
  1,18,-72.2349,42.555789,ATHOL,15 2,21,-72.2349,42.582768,ATHOL,15
  2,23,-72.2349,42.600754,ATHOL,15 2,24,-72.2349,42.609747,ATHOL,15
  2,34,-72.22273,42.582768,ATHOL,15 2,35,-72.22273,42.591761,ATHOL,15
  2,36,-72.22273,42.600754,ATHOL,15 2,37,-72.22273,42.609747,ATHOL,15

I am graphing column 3 as the X-coordinate and column 4 as the Y-coordinate.

However, I would really like to color-code these by column 1 so it is easier to see the groupings on my map. For instance, any coordinate where column 1 = 1 can display a random color (ex. green) and coordinates where column 1 = 2 can display another random color (ex. red).

I kind of found a way to do this, but it is extremely tedious as I have a couple thousand data points and a hundred or so items in column 1.

What I have been doing is going to Properties > Style > changing from Single Symbol to Categorized > setting column to 1 > setting a random color ramp > then add categories and setting their value manually from 1, 2, etc.

Is there an automated way to do this, so I don't have to set hundreds of categories and enter values manually?


texts:
I have the below code that on a short story calculates Leaf Area Index(based on work of Schiffman et al., 2008). When i want to export the LAI collection I transform it in a list to export all the elements of the list. I used this method  to export collection images in the past and it work as smooth as possible everytime.

The below code is for running the LAI script:
https://code.earthengine.google.com/?scriptPath=users%2Fgeorgeboldeanu20%2Fscripturi_bune%3ALAI_L8

Right now I have problems with the exports, it exports it into tiles
(2-3 maybe 4 tiles), it exports parts of the whole images(some are missing), it exports tens of small images that are emptyall fixed in size. Someone's been dealing with this issue lately??

I didnt find any answer for this, but i saw some google earth engine apps down, maybe there are some cloud issues?


texts:
Trying to do some research into drought patterns of the US using KMZ files from droughtmonitor.unl.edu.

My problem is that there are too many files to work with. 857 total.
Records are every 7 days from 01/04/2000 - 03/22/16

I would like to merge each year into one KML. Doing this in Google Earth from importing all the separate files is very time consuming and causes GE to crash about halfway through. Also I would like to clean up some folder names. Each folder is 66 characters and is repeated in the subfolder.

How would I merge multiple KMZ/L files into a single KMZ/L file and clean up the folder names?

Example:

From

Top folder: U.S. Drought Monitor Released: 2000/01/06 Valid 7am EST 2000/01/04
  Subfolder: U.S. Drought Monitor Released: 2000/01/06 Valid 7am EST 2000/01/04
Top folder: U.S. Drought Monitor Released: 2000/01/06 Valid 7am EST 2000/01/11
  Subfolder: U.S. Drought Monitor Released: 2000/01/06 Valid 7am EST 2000/01/11

To

Top Folder: Merged 2000 Drought Maps
  Subfolder: 2000/01/04
  Subfolder: 2000/01/11

Here is a link to a drought map KMZ.
http://droughtmonitor.unl.edu/data/kmz/usdm_20000104.kmz


texts:
I've just installed QGIS 2.8 on the new Fedora 22 (64) that came out some days ago. Awkwardly, I can't find the GDAL Georeferencer that, to my knowledge, is supposed to come with every installation. Searching for it in the  Plugin Manager doesn't give any results so I'm unable to activate it thereby adding it to the Raster menu. Re-installing it all is no good either.
Any idea what might be the cause for this? Just two days ago, everything worked perfectly fine (2.8 running on Fedora 21). Or is it just me being lost in the woods?

Recently, I downgraded back to Fedora 21 (less buggy at the moment) but am now facing the exact same problem. Is the Georeferencer included somewhere else by now?


texts: Methods to accurately georeference geographic map?

texts:
I'm writing an application that allows the user to sample points and perform georeferencing

In order to test it I have a large scanned map (16k x 13k) pixels covering the area 31 north to 34 north and 34 east to 36 east with a geographic grid

I sampled the 4 corners of the map + 2 points on the diagonal
and applied an affine transformation :

[x']   [A B C]   [x]
[y'] = [D E F] * [y]
[1]    [0 0 1]   [1]

using the 6 points, and least squares I managed to find [A B C D E F], but applying it gave me an inaccurate transformation...

The Longitude is very accurate all along the picture (> 1 arc second), but the latitude deviates sometimes up to 15 arc seconds (~ 400-500 meters)

I assume that no matter how many points I will sample the affine transformation won't cut
it...

What is the appropriate transformation here?

Trying to project the geographic coordinates

I projected the sampled points to UTM first, and then calculated the affine transform (again 6 points)

The results are pretty much the same - Accurate near the control points but, drops significantly when I get farther

Using 2nd Polynomial Transformation

I was able to get better results by using a 2nd order polynomial transform instead of an affine one, but I still get deviation of 4-5 arc seconds

The transformation I used -

[x']  =  [A B C D E F] * [x*x]
[y']     [G H I K L M]   [y*y]
                         [x*y]
                         [ x ]
                         [ y ]
                         [ 1 ]

Where A-M are the transformation parameters {x',y'} are the geographic coordinates and {x,y}
are the screen pixel coordinates.

since there are 12 parameters it requires 6 points to be sampled...

I suppose I could try to do a 3rd order transform (20 parameters - 10 control points needed)
or even higher but I wonder if there isn't a better way?


texts: error "pointer hDs is NULL in GDAL GetProjectionRef" when using gdal_rasterize

texts:
I want to convert a .shp to .tif.
So I find the tool gdal_rasterize may help.
But when I input the file it always meet this error:

pointer hDs is NULL in GDAL GetProjectionRef

I have tried different versions of GDAL(1.11.0 and 2.1.0),same error.
The test data is from http://qgis.org/downloads/data/qgis_sample_data.zip

1)Is this wrong because of the data?or the wrong syntax?

2)What is the correct syntax for using gdal_rasterize()?

3)Is there a better way to convert shp to png(or tif)?coding by myself is ok but I use Java.


texts:
The title says it all. I'm using GeoExt as part of OpenGeo Suite. I would like to display the current mouse coordinates.

Update: I was getting confused as to what parts of the code are GeoExt and which parts are OpenLayers.

Here's some code. mapPanel is a GeoExt object. mapPanel.map is an OpenLayers object.

var mapPanel = new GeoExt.MapPanel({
    // etc.
});

mapPanel.map.addControl(new OpenLayers.Control.MousePosition());


texts: Insert n points equally distributed into many polygons

texts:
I know the "Random Points inside Polygons" tool in QGIS. It generates the below output.

What I want to achieve however are 8 equally distributed points for each of the polygons. I thought of dividing the polygons Voronoi-style by the random points and then generating new (equally distributed) points, but I don't seem to be able to do so. Using "Voronoi Polygons" in the toolbox, I only get a Voronoi diagram that's bounded by a rectangle. No way to divide the whole polygon by it:

How can I get 8 equally distributed points inside each of the polygons?


texts: Where can I find a good Sextante Modeller tutorial?

texts: I'd like to enrich my knowledge of the Sextante Modeller. Do you have any good suggestion of a book/workshop etc. to dig into the modeller?


texts:
I'm trying to create a shapefile of a state. I have a map of US zipcodes that I downloaded from the US Census and I want do create a shapefile of zipcodes in Florida.

I ran this query which produces a shapefile that I expected, however when I examine the attribute table in QGIS, I see that there are duplicates of the zip code. The same thing occurred when I modify the query to intersect a specific state using a shapefile of US DMAs (Flordia has 11 DMAs, but the attribute table has 1,107 rows, mostly duplicates). Answers regarding how to fix duplicate rows in an attribute table has come up unsatisfactory. An example of my output is here:

My query is below:

"SELECT geoid10, geom FROM (
    SELECT uszips.geoid10, (ST_Dump(ST_Intersection(uszips.geom, fl.geom))).geom
    FROM uszips, fl
    WHERE ST_Intersects(uszips.geom, fl.geom)) t
WHERE ST_GeometryType(geom) = 'ST_Polygon';


texts: Geoserver data directory on a network drive?

texts:
is it possible to set the data directory of geoserver on a network drive? The data directory path of my existing installation is the default location (C:\Program Files\GeoServer 2.6.1\data_dir). I copied the folder on a network drive (a mapped drive: I:\GeoServer_Data_Dir (read and write permission for the used Windows-User)) and change the environment variable (GEOSERVER_DATA_DIR) to the new location. After restarting the geoserver service he didn´t use the "new" data dir. Instead he uses C:\Program Files\GeoServer 2.6.1\webapps\geoserver\data as the data directory. In addition i installed GeoServer 2.8.3 on this Computer and during the installing i set the data directory path to the existing folder on the network drive. It said that the path contains a valid data directory but the GeoServer 2.8.3 also uses the webapps\geoserver\data-folder in his home-directory.

Has anybody an idea how to fix this problem?


texts: Project with multiple MapSeries - script to select specific series and range

texts:
Currently working a project that contains 9 map series with over 22,000 pages total - split across datum and scale. Want to script a solution that allows us to export pages as individual PDFs by selecting a specific series and range - ultimately to allow us to have multiple PC's exporting maps at once.

All the examples on the ArcGIS Pro help appear to be referencing projects with a single Map Series so I'm a bit stuck how to approach. So far I've written a small script to check for the presence of the series and return the page count:

    print("Layouts:")
    for lyt in aprx.listLayouts():
    print(lyt.name)
    ms = lyt.mapSeries
    if ms is None:
            print('\t'+"None available")
    else:
        if ms is not None:
            print('\t'+"Page count:",ms.pageCount)

I can export a PDF based off a selection on the index layer - however the problem I am finding is I don't know which map series the selection is operating on. Anyone have to suggestions to specify which map series? Would ideally prefer to go with ranges than selecting attributes (ie 1 - 500, 501 - 1000 etc)

#make a selection based off index feature
if not lyt.mapSeries is None:
    ms = lyt.mapSeries
    if ms.enabled:
        ms = lyt.mapSeries
        indexLyr = ms.indexLayer
        arcpy.SelectLayerByAttribute_management(indexLyr, "NEW_SELECTION", "ZONE = '56'")
        ms.exportToPDF(Output + "\\TEST.pdf", "SELECTED")


texts: Is it possible to write a server to handle OpenLayers.Protocol.HTTP({url: "/cgi-bin/fs/workspace/" ...}) from FeatureServer

texts:
I am a newbie to geoserver and openlayers. Recently, I am stuck in saving changes(drawing polygons and points) of the feature vectors to a local(geoserver/www/...) .json file.

The tutorial I followed is here.
But it said:

  "Of course you also need a program on your server that knows what to
  do with the transactions sent by the protocol/strategy. The server
  side has nothing to do with OL, but you can't demonstrate
  functionality without it, so I'll use a workspace I play around with
  (which uses FeatureServer, which sadly has not been much supported in
  recent years). This is only for small-scale test purposes, and any
  features entered in it may disappear at any time."

Besides, http://featureserver.org/ cannot be visited. So I wonder if it is possible to write a server program myself and replace the url with "http://localhost:XXXX/MyWebServer/index.jsp" under the OpenLayers.Protocol.HTTP. Then the index.jsp file will handle the GET/PUT request and response to send or save the json file.

Since I am a beginner, correct me if anything I mentioned is wrong. That would be appreciated.


texts:
I need to create my own mapsforge map file...

My system is Ubuntu 14.04.1 LTS
I've installed osmosis via apt
Since link in official documentation is not working i've found plugin file here and downloaded it
Following this documentation I placed file in ~/.openstreetmap/osmosis/plugins/
Downloaded osm.pbf file from geofabrik
Run

osmosis --read-xml file=poland-latest.osm.pbf --mapfile-writer file=test.map

And i' getting this error:

lis 17, 2014 10:40:02 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Osmosis Version 0.40.1
lis 17, 2014 10:40:02 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Preparing pipeline.
lis 17, 2014 10:40:02 AM org.openstreetmap.osmosis.core.Osmosis main
SEVERE: Execution aborted. org.openstreetmap.osmosis.core.OsmosisRuntimeException: Task type mapfile-writer doesn't exist.
        at org.openstreetmap.osmosis.core.pipeline.common.TaskManagerFactoryRegister.getInstance(TaskManagerFactoryRegister.java:60)
        at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.buildTasks(Pipeline.java:50)
        at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.prepare(Pipeline.java:112)
        at org.openstreetmap.osmosis.core.Osmosis.run(Osmosis.java:86)
        at org.openstreetmap.osmosis.core.Osmosis.main(Osmosis.java:37)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:483)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launchStandard(Launcher.java:328)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launch(Launcher.java:238)
        at org.codehaus.plexus.classworlds.launcher.Launcher.mainWithExitCode(Launcher.java:408)
        at org.codehaus.plexus.classworlds.launcher.Launcher.main(Launcher.java:351)
        at org.codehaus.classworlds.Launcher.main(Launcher.java:31)

Seems like osmosis don't see the plugin, but i've done everything according to docummentation... I've also tried to place plugin file in other localization (/usr/share/osmosis, /usr/share/osmosis/plugins), but still not works..

What am I doing wrong?


texts: Converting raster (array) to NetworkX graph and do * search?

texts:
I want to use a slope-raster for a a_star path analysis in NetworkX.
I am using Python for this project.

My test raster looks like this

First I read in the raster with GDAL

input_raster = "slope5.tif"
raster = gdal.Open(input_raster)

Next I read the raster as an array

bandraster = raster.GetRasterBand(1)
arr = bandraster.ReadAsArray()

The array looks like this

>>> arr
array([[ 15.26072407,  16.01073837,  27.81685448,  24.49266815,
         32.94010162],
       [ 26.0834446 ,  21.01086426,  30.02531815,  25.44641304,
         19.48436928],
       [ 21.98579979,  19.06574059,  24.44838905,  26.67520905,
         26.61030769],
       [ 20.42549133,  28.58900261,  17.00629425,  31.25236893,
         32.34066772],
       [  8.33120537,  24.303339  ,   7.85025167,  35.50911713,
         29.56856346]], dtype=float32)

Then I create from the array a graph

G = nx.DiGraph(arr)

As far as I understand it my graph looks like this

(0,0) (0,1) (0,2) (0,3) (0,4)
(1,0) (1,1) (1,2) (1,3) (1,4)
(2,0) (2,1) (2,2) (2,3) (2,4)
(3,0) (3,1) (3,2) (3,3) (3,4)
(4,0) (4,1) (4,2) (4,3) (4,4)

So for example the weight of G[3][2] is {'weight': 17.00629425048828}

How do I now find the A star path?

I tried the following to find the way from (1,0) to (3,3)?

nx.astar_path(G,(1,0),(3,3))

But that only gave

Traceback (most recent call last):
  File "<pyshell#44>", line 1, in <module>
    nx.astar_path(G,(1,4),(4, 2))
  File "C:\Python27\lib\site-packages\networkx\algorithms\shortest_paths\astar.py", line 108, in astar_path
    for neighbor, w in G[curnode].items():
  File "C:\Python27\lib\site-packages\networkx\classes\graph.py", line 319, in __getitem__
    return self.adj[n]
KeyError: (1, 4)

My nodes and edges look likes this:

>>> G.nodes()
[0, 1, 2, 3, 4]
>>> G.edges()
[(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (1, 0), (1, 1), (1, 2), (1, 3), (1, 4), (2, 0), (2, 1), (2, 2), (2, 3), (2, 4), (3, 0), (3, 1), (3, 2), (3, 3), (3, 4), (4, 0), (4, 1), (4, 2), (4, 3), (4, 4)]

Is this the right approach?
How do I find the path?
Did I mix up something with nodes and edges?


texts: Merging address fields, inserting ", " where needed

texts:
I have a large address database point FC in FGDB format. Each address record is split into multiple fields like 'organisation', 'building name', 'building number', 'post town', 'post code' etc. Not every field is populated with data, there are many null values (not every address has an 'organisation' name for example).

I've created a new text field called 'Full_Addr' which will merge all the address fields into one. So far I used this VB script in field calculator:

Full_Addr=
[OrganisationName]& ", " & [BuildingName]& ", " & [BuildingNumber]& ", " & [Thoroughfare]& ", " & [PostTown]& ", " & [Postcode]

However this results in entries like:
, , 1, Example Street, Example Town, AB1 2CD
or
Example Organisation, , , Example Street, Example Town, AB1 2CD

Is there a smarter, more efficient way of creating a merged address field without the double, triple or quadruple ", "?


texts: QGIS 2.8 Convert a vector to a raster layer

texts:
When I change a vector layer to raster using raster button->converse->Vector rasterization in QGIS2.8. I get only a black picture.
Why is that so and how can I convert a vector layer to a raster layer.