texts: Updating xy fields in attribute table automatically when moving points to new locations in ArcMap?

texts:
I work on ArcGIS Desktop 10.3.1 and try to move several points so they are corresponding simultaneously to the coordinates value in the attribute table immediately.
I know this question repeats in this site and I read those questions and their answers:

Update xy fields when I move points?
Move points to XY coordinate
Dynamically update points with xy coordinates in table
Updating Location of Point Features using XY Fields?
https://community.esri.com/groups/technical-support/blog/2011/12/06/how-to-update-the-location-of-a-point-feature-and-its-xy-fields

I seek for a solution with no python\any code or using the calculate geometry

I also tried to create a feature class in GDB from the shapefile and in order to get automatic update XY fields in the attribute table

Is it possible to do it in some way with GDB or with another way?


texts: Why can't I open QGIS projects via windows explorer?

texts:
We have upgraded from QGis 1.8 to QGis 2.2. Qgis 2.2 is first installed, then 1.8 removed. Since then it is impossible to open .qgs-projects directly from windows explorer.
We have checked explorer, file extensions and relinked to C:\Program Files\QGIS Valmiera\bin\qgis.bat . But it doesn't help.

We get the error 'qgis_core.dll cannot be found, please reinstall'. But that is rather drastic.

I use windows xp.


texts: Merging multiple GeoJSON files into one file with layers

texts:
What is the best way to merge several GeoJSON files of different types into a single GeoJSON but where each original GeoJSON is its own data layer within the final output?
Been playing with QGIS and can't figure out a way to get this result without combining all the data together into one layer (which when I have polygons and points I can not do).
Can most likely write a script to read in and convert all the JSONs and then write out to a single file, but just wondering if there is a better/built in way of doing this.


texts:
I'm using ArcGIS 10.2 Desktop on a Windows platform.

I have georeferenced about a dozen USGS topographic raster maps from around 1900.  Each georeferenced map is an now an .mxd file.  To the best of my limited knowledge there are no layers associated with any of these maps.

I would like to tile these into a single map [to which multiple layers would be added] but cannot find a way to do this.  Perhaps I need to start over using a different technique?

Would be VERY grateful for any suggestions on how to solve this problem.


texts: What is the significance of the SDE lineage length?

texts:
ORA-55378 left a great answer to my question:

Does a high number of SDE_state_lineages slow performance?

But, I realized I don't really understand what lineage length means and how it affects performance? Using the example on this ESRI technical article:

VERSION NAME              LINEAGE LENGTH
------------------------- --------------
TOM.WORKORDER_I0C6531     21
TOM.EDITOR_TB             25
ERIK.WORKORDER_I910313    78
BRENT.DESIGN_F831         92
BRENT.DESIGN_F831_ALT_1   94
BRENT.DESIGN_F831_ALT_2   95
SDE.MAPPING_DEFAULT       138
SDE.LANDBASE              169
SDE.DEFAULT               249
ERIK.EDITOR_EH            275
STEVE.EDITOR_SM           289

What do these lineage lengths actually mean? What is the significance of a length of 289?


texts:
I'm looking for software that can perform local adjustments in DEMs with tools similar to paintbrushes used in graphics editing programs but streamlined to respect geomorphology.

To give you an idea of what I have in mind check out this video of SimCity 4 terraforming mode. I would give a lot to have such a tool in ArcGIS environment.

Any local DEM modification in ArcMAP is a hassle.

Let's say I have a DEM with river bathymetry imprinted in it but there is a river port along that wasn't covered and has the height of the water surface (that's the usual case with LIDAR derived DEMs - water bodies are filled).

The only information about the harbour is the maximal depth. So to transfer this information to my DEM, I erase the entire harbour area (null values), imprint a small raster with maximal depth in the center of the harbour basin and then I fill the remaing null cells with a moving window filter.

The results of this operation are usually sufficient but the whole process is complicated. If only I had a "lower terrain paintbrush" :)


texts: How to create custom discrete color scheme for polygon attributes in ggplot2?

texts:
Back again with more questions. :)

Data: Zipped ea.ml1 shapefile can be found here, and the country outline can be downloaded here at the GADM site. I've been using the Ethiopia ESRI shapefile.

As per suggestions on a previous question of mine, I have decided to try to transform my ongoing map project/test to ggplot2 for ease of use. You can check out my full original non-ggplot code here if you're into that kind of thing. I have done a tiny bit of mapping with ggplot2 in the past, so I'm kinda-sorta familiar with it, but I've certainly not mastered it. I'm wondering how to change/customize my color scheme for the ML1 attributes of my Ethiopia polygon. It should look like this (which I did with base plot functions, as seen in my linked code):

Using this brief tutorial from GitHub, I managed to finagle with my code a bit to create this:

library(ggplot2)
library(maptools)
library(rgdal)
library(rgeos)
library(raster)
library(plyr)

setwd("D:/Mapping-R/Ethiopia")

ea.ml1 <- readOGR(dsn = "D:/Mapping-R/Ethiopia", layer = "EA201507_ML1")
eth <- readOGR(dsn = "D:/Mapping-R/Ethiopia", layer = "ETH_adm0")

#remove all unwanted columns from eth shapefile
eth <- eth[, -(2:67)]

#define color scheme for ml1 & ml2 features
ea.ml1@data$COLOR[ea.ml1@data$ML1 == 1] <- "gray"
ea.ml1@data$COLOR[ea.ml1@data$ML1 == 2] <- "yellow"
ea.ml1@data$COLOR[ea.ml1@data$ML1 == 3] <- "orange"
ea.ml1@data$COLOR[ea.ml1@data$ML1 == 4] <- "red"

#clip ethiopia to ml1 background
eth.clip1 <- intersect(ea.ml1, eth)

#identifies attribute rows, converts to data frame, joins points to attributes
eth.clip1@data$id <- rownames(eth.clip1@data)
eth.f <- fortify(eth.clip1, region = "id")
eth.df <- join(eth.f, eth.clip1@data, by = "id")

ggplot(eth.df) +
  aes(long, lat, group = group, fill = ML1) +
  geom_polygon() +
  geom_path(color="white") +
  coord_equal()

which resulted in this:

How do I customize my ggplot2 color scheme (which is actually defined in the code, not sure if that's useful or not) to match the first map image that I posted (where an ML value of 1 is gray, 2 is yellow, etc.)? The linked tutorial simply lists its color scheme as "Utah Ecoregion" with no details about where it came from/what it looks like. I also tried the solution listed here to no avail; R just spit nonsensical errors at me.

 ggplot(eth.df) +
  aes(long, lat, group = group, fill = ML1) +
  geom_polygon() +
  scale_fill_manual(values = (1 = "gray", 2 = "yellow", 3 = "orange", 4 = "red")) +
  geom_path(color="white") +
  coord_equal()

gives me these errors:

  Error: unexpected ',' in: "  geom_polygon() + scale_fill_manual(values = (1 = "gray","

  Error in geom_path(color = "white") + coord_equal() : non-numeric argument to binary operator`

So, how do I customize the color scheme without screwing everything else up?


texts: Keep "no data" tiles with r.reclass rules

texts:
Working in QGIS 2.14.0-Wien.

I want to reclassify a raster with the r.reclass GRASS 7 tool, but want to keep the "no data" tiles as "no data".

I understand that I need to create a text "rules" file to inform r.reclass of what classes to create with lines such as "XXXX thru XXX = 1", but am not sure of the syntax to preserve the "no data" tiles.

How do I do that?


texts: Filling voids in DEM SRTM 1-arc second using QGIS by resampling and inserting another DEM

texts:
I'd like to use the DEM SRTM 1-arc second (Version 4, 30m) from USGS for Europe. Unfortunately there are some voids. My idea is to fill those voids using the void-filled version 3 (90m) by resampling it to 30m and inserting the missing data in the version 4 DEM.

As I'm new to QGIS I don't know the proper steps. I already got a 30m map by reprojecting the dimensions times three but I'm not sure if that's the right thing to do and I have no idea how to fill the voids afterwards.


texts: arcpy.da.UpdateCursor - how to make multiplication of each column with the selected column, for each row in particular?

texts:
If I have a table with a number of columns [col1, col2, col3 ... col_n] how to perform        multiplication of each column with one selected column, for example:

multiplication = col1 * col_def, col2 * col_def, col3 * col_def ... col_n * col_def

and then to determine the sum, for example:

sum_multiplication = Sum (multiplication)

I try as follows, but does not go:

with arcpy.da.UpdateCursor(fc, (fieldNameList)) as cursor:

   for row in cursor:

      for i in range(len(fieldNameList)):
         hi = sum(row[i] * row[14])
         row[15] = hi
      cursor.updateRow(row)

Does anyone have an idea how to solve this problem?


texts: Geodesic measurements for short distances throughout US?

texts:
My ArcGIS experience has been limited to working on large-scale maps with planar projections. I'm moving to a project in which I want to calculate the distance between pro sports facilities throughout the U.S. and their city's respective downtown/central business district (both defined as geocoded points).

My primary interest is in relative distances, and my question is whether using a geographic coordinate system for US map points and using "Near" to calculate geodesic distances between stadia and downtowns is a sound approach to avoiding issues that would arise from calculating distance based on a projection of the US?


texts: Changing data type in ArcGIS attribute table?

texts:
How can I change the data type in an attribute table in ArcGIS for multiple fields?

For example, if the fields have been exported as a string from Excel, and I would like to use them as numerical values?


texts: Creating polygons with fixed extent based on point which represents S/W corner using QGIS

texts:
I have a point data set which contains the information of workplaces per hectare.
Each point represents the South-Western corner of the surveyed hectare. Not the entire map is covered with points, only there where people are working does a point even exist.
I want to create now a polygon shapefile with many rectangular polygons (100*100m) and the information from the point.
In the image you can see how the point relates to the desired polygons. The ones scribbled in black are not necessary, as there are no workspaces.
I am using QGIS.


texts: Editing WFS-T layers from ArcGIS Server in QGIS?

texts:
I'm trying to edit layers from a service in ArcGIS Server through WFS-T in QGIS. But the AGS only allows writing in version 1.1 of WFS, while QGIS only connects to the WFS Server in version 1.0, making it impossible to start editing the layer.

Does anyone know any way to solve this issue, or the development of QGIS walks into a future update of its WFS module to version 1.1?


texts:
Is there any description of OGC WMS specification in OpenAPI / Swagger format ?

I found an initiative for WFS protocol, but nothing on WMS.


texts: Adjusting pixel values of a raster based on geopandas colums or fields

texts:
I have a geodataframe with geometry column containing Polygons and 5 fields or columns('Crop1', 'Crop2', 'Crop3', 'Crop4' and 'Crop5') as well as a raster image. There are 159 rows in the geodataframe which also means there are 159 Polygons(in the geometry column). I am supposed to "create a masked version of the raster I have where pixels falling within one of the fields are set to 1 and pixels outside the fields are set to 0". I thought about using rasterio.plot.adjust_band, but after reading it over and over, I don't think that is the solution. The geodataframe and the raster image have the same CRS.
Any ideas?


texts: Seeking tool in QGIS similar to ArcGIS Combine?

texts:
I am wondering if there is a way/tool in QGIS (2.6) could produce similar output as ArcGIS "Combine" tool does. I am trying to integrate values from overlapping rasters, and then create a single raster with unique value for each unique combination of input values. If there is a tool in QGIS could be used to get that will be the best. If not, should I using a script to automate reading through each pixel values? I am afraid it will take a large amount of time to processing national-wide data.

Here is the link to document for the "Combine" tool:
http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/combine.htm


texts:
I am quite new in QGIS and I really need your help.

I have a map of Kenya divided into 68 districts and I have coordinates of 30.000 schools. The points have one variable besides the latitude and longitude data called "distance to nearest school".

First I imported the map. Then I imported the spreedsheet and transformed it to a SHP file like the map. So now I have two SHP files.

The points have three variables: "Latitude", "Longitude" and "Distance to nearest school". I want to have a fourth variable called "District" and take the mean of the "Distance to nearest school" for each district.


texts: Adding another style to an OpenStreetMap server using Mapnik, PostGIS, etc

texts:
I have created an OpenStreetMap server using the guide Manually building a tile server (12.04). It's very good and works like a charm. I have also been able to add a new theme using TileMill and osm-bright.

However, I now want to be able to switch between two themes, osm-bright and my own osm-night. Is this an Apache question or is there a way to modify Mapnik to be able to get a parameter from, for example, openlayers to switch between these two themes?


texts: GSEA: unable to find an inherited method for function ‘keytypes’ for signature ‘"standardGeneric"’

texts:
I can't understand this problem.  Even when I use the data and code offered at https://learn.gencore.bio.nyu.edu/rna-seq-analysis/gene-set-enrichment-analysis/ , the problem also existed.

I did the analysis follow the official illustrate. As bellows:
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(org.Dm.eg.db)
library(DOSE)
library(biomaRt)

df = read.csv("~/lineagepaper/drosphila_example_de.csv", header=TRUE)

original_gene_list <- df$log2FoldChange

names(original_gene_list) <- df$X

gene_list<-na.omit(original_gene_list)

gene_list = sort(gene_list, decreasing = TRUE)
head(gene_list)

gse <- gseGO(geneList=gene_list,
             ont ="ALL",
             keyType = "ENSEMBL",
             nPerm = 10000,
             minGSSize = 3,
             maxGSSize = 800,
             pvalueCutoff = 0.05,
             verbose = TRUE,
             OrgDb = organism,
             pAdjustMethod = "none")


texts: QGIS: Zoom to selected feature via plugin doesn't works in first attempt

texts:
I recently migrated from QGIS 2.18 to QGIS 3.6. In the plugin, I load a vector layer and then select and zoom to a particular feature.
This worked fine in v2.18 but now on 3.6, what happens is the layer is loaded, feature is selected but zoom to that feature doesn't work in first attempt.

But when I open my plugin again in the same session (without closing the previous instance), it works fine.

Code used to zoom to feaure:

def zoom_to_tile(self):
    print("zoom_to_tile")
    print('zooming to tile:', self.current_tile)
    box = QgsGeometry.fromWkt(str(self.current_tile[2])).boundingBox()
    box.grow(.002)
    # self.iface.mapCanvas().setExtent(box)
    # self.iface.mapCanvas().zoomToSelected()
    # self.iface.mapCanvas().refresh()
    print("selected ",self.mortonTileVLayer.selectedFeatures()[0][0],
          " iface:",self.iface.activeLayer().selectedFeatures()[0][0])
    canvas = self.iface.mapCanvas()

    canvas.zoomToSelected(self.mortonTileVLayer)

    canvas.refresh();
    print("zoomed")

I need help to get it working first time around.


texts: Plotting a half degree grid on a map and extracting species information using R

texts:
I have a dataset which looks like this(subset):

Latitude        Longitude   Species
-34,68871595    20,0331831  Ischnura senegalensis
-34,59082383    20,05421162 Ceriagrion glabrum
-34,59082383    20,05421162 Crocothemis erythraea
-34,59082383    20,05421162 Ischnura senegalensis
-34,63333333    20,05       Africallagma glaucum
-34,53333333    20,05       Africallagma glaucum
-34,53333333    20,05       Ceratogomphus pictus
-34,53333333    20,05       Crocothemis sanguinolenta
-34,30          18,45       Orthetrum julia
-34,35          18,45       Orthetrum julia

First I plot the data on a map of Southern Africa. What I want to obtain as output is a dataframe where for each cell of my grid I obtain the number of individuals within the cell, as well as number of species. The cell size I am wanting is a half degree cell which is align with the centre of the latitude and longitude of my map i.e one cell centre will be on -34(lat) and 20(long)

My code insofar-
###plotting data points
lonlat<-cbind(Longitude,Latitude)
library(sp)
pts<-SpatialPoints(lonlat)
showDefault(pts)
crdref<-CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
pts<-SpatialPoints(lonlat,proj4string = crdref)
Sppnames<-SA.ODA$Species
df<-data.frame(lonlat,Sppnames)
ptsdf<-SpatialPointsDataFrame(pts,data = df)

###Plotting a Africa base map
mapofSA<-readRDS('/Users/Ashleigh/Desktop/gadm36_ZAF_0_sp.rds')
proj4string(mapofSA)<-CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')


texts: ST_Union or ST_Collect on PostgreSQL union result?

texts:
Newbie to PostGIS/Postgres. I have the follow query which returns me 2 GEOMETRYCOLLECTIONS. Is there a way within my current query to return one polygon?

SELECT ST_AsText(ST_Collect(ST_Simplify(geom,.005))) AS singlegeom
FROM public.mzones
WHERE id = ALL(ARRAY['GMZ730'])
UNION
SELECT ST_AsText(ST_COLLECT(ST_Simplify(geom,.005))) AS singlegeom
FROM counties
WHERE fips IN (01031,01045,01061,01067,01069,12005,12013,12037,12039,12045,12059,12063,12073,12077,12129,12131,12133,13061,13087,13099,13201,13253)
;


texts: How does Mbtiles store raw image data in the tile_data field?

texts:
Looking at the MBtiles spec, the tile_data field of the tiles table is just a blob of a tile image.

With this information, I expected to be able to take a data URI, string it of everything except the base64 encoded string, and decode it for storage in the tile_data field.

For instance:

var dataURI = "data:image/jpeg;base64,[some_encoded_data]";
var unencodedBytes;
if (dataURI.split(',')[0].indexOf('base64') >= 0) {
    unencodedBytes = atob(image.split(',')[1]);
} else {
    unencodedBytes = unescape(dataURI.split(',')[1]);
}

At this point, my unencodedBytes variable should have raw jpeg image data, and as such, I should be able to insert it into a tile_data field, correct?  Is there an additional step I need to take?  The data looks like belongs there, based off my observations of other valid mbtiles files, however all the mbtiles viewers I try to view this either do not show the data, or say they cannot figure out what type of data it is (jpeg or png).

What am I doing wrong?  This is all being run in javascript in a browser, and I am using WebSQL to create the mbtiles file, since it is (supposedly) compatible.

Thanks in advance


texts: Count all records X distance from a given point using SQL in CartoDB

texts:
I'm trying to use CartoDB to Count all of the records a given distance from a selected location.

I have the lat/long for the selected location (call this a store). Now I need to find all of the records (representing customers) that are a given distance (say, 8046 meters/5 miles) from the selected location.

What SQL query should I be using? I've tried some basic things like this, but they didn't work:

'SELECT COUNT() from my_table_name WHERE the_geom && ST_SetSRID(ST_Point(39.8282, -98.5795), 8046)'

or this:

SELECT COUNT() from my_table_name WHERE ST_INTERSECTS(the_geom,CIRCLE(LATLNG(39.8282, -98.5795),8046))

What query do I actually need? All I want to do is count the records that meet this criteria. Step 2 would be to filter the points on the map to show only those that meet the criteria, but one thing at a time.


texts: "Please ensure that the GDAL_DATA environment variable is set and try again" on Kubuntu 12.10

texts:
I'm having a problem loading the 'convert os gml to ogr' plugin in QGIS 1.8 running on (k)ubuntu 12.10.

Basically, clicking the plugin prompts the error:

"Please ensure that the GDAL_DATA environment variable is set and try again"

I've tried setting the environment variable via instructions here: (http://trac.osgeo.org/gdal/wiki/FAQInstallationAndBuilding#WhatisGDAL_DATAenvironmentvariable) with the terminal command "export GDAL_DATA=/usr/share/gdal/1.9"
but to no avail (i'm not even sure this is the right course of action...)
Any suggestions?


texts: Reproject shapefile in QGIS using "save as" return weird result

texts: I need to reproject a shapefile from a geographic cooridante system (ESPG: 4326) to a projected coordinate system (EPSG: 54004). Obviously something went wrong but I'm not sure what.


texts: Change directory where Layer Preview generates shapefiles in GeoServer

texts:
How to change directory where GeoServer puts temp shapefiles generated for export via Layer Preview?

Currently it is putting them in .../tomcat/temp/XXX and I would like to put them to another directory lets say ./usr/myTempDir?


texts: Geodjango, Openlayers 2: Override default fillColor in OpenLayers.Feature.Vector.style

texts:
Update

I would certainly hard-code the source file
and change it as follows:

OpenLayers.Feature.Vector.style = {
    'default': {
        fillColor: "blue",
        ...
    },
    ...

... but this is rather a hack I would like to a avoid.

I am trying to change the fillcolor of a point shown in a django admin map.

I have already managed to change the color of a new point when I click on the map, following this solution, but I am unable to understand how I could change the color of an existing point that programmatically is inserted into the map when fetched from the database.

I have found several posts that describe how to override the StyleMap of a vector layer or how to dynamically style a map, but I do not understand how I could change the color of the dot when it is already on the map.

I attach a picture. The yellow dot is already in the database.

How can I paint it blue?


texts: How to hierarchically query PostgreSQL data coming from OpenStreetMap

texts:
I have loaded OpenStreetMap data into a PostgreSQL database using Osmosis. This is done by using the snapshot schema.

I filtered the data and loaded only: country & city information

osmosis --read-pbf planet-latest.osm.pbf --tf accept-nodes place=country,city --tf reject-relations --tf reject-ways --lp --write-pgsql host="localhost" database="testing" user="postgres" password="s5mLzcqy"

This data is correctly loaded, but now I want to query this data. The output of this query should be the name of the city and next to it the name of the country.

I already have this query which extracts the name of the cities:

SELECT
    tags->'name' AS name
FROM
    nodes
WHERE
    tags->'place' IN ('city')
;

Can someone tell me how I can also get the names of the country the city is in?

Update:
I tried to work with the ST_Contains() function, but I can't get it to work. Is there someone that has experience in working with this function?


texts: Working out the area of certain class types within a polygon

texts:
I am trying to use a DEM to determine land management units based on combinations of slope and aspect.

As per the task I have created three individual boolean 'Class' rasters which have specific aspect/slope details. I also have a shapefile of land regions - which has area Ha data contained per region.

I want to work out the area of each "Class" within each region. I know this shouldn't be that difficult but I can't work out which tool I should be using.

Thanks,


texts: How to fix ERROR 6: Unable to load PROJ.4 library (proj.dll), creation of OGRCoordinateTransformation failed issue

texts:
I have installed GDAL in my system. I tried a simple program which uses spatial reference functions. When I tried to execute that program i am facing "ERROR 6: Unable to load PROJ.4 library (proj.dll), creation of
OGRCoordinateTransformation failed" issue. I have installed proj.4 also. I have also added proj.4 path in environment variable. Still I am facing this issue. How to fix this?


texts: select which is having specific double value

texts:
Title: Trying to Select which is having double space

I have data like below:

field name - Address:

Winding Springs Dr
Winding Hollow Dr
Windsong Ln

in above data i have to select double space values those are like 1 & 2

i'm following:

Adding dummy field - coping data into new field
Removing spaces in new field
Then finding string lengths of new and old field

finally which is having two or more length differ then that is my output

It is too lengthy a process so can you advise a better solution?


texts: How to prevent zoom-in when clicking in the Leaflet table of contents window?

texts:
The Leaflet grouped layer control plugin has a demo page, where double-clicking within the Table of Contents window has no effect on the map. This sample uses Leaflet 0.73.

When using this plugin with Leaflet 1.x, I'm finding that when I double-click within the table of contents window the map zooms in. You can test this in the simple JSFiddle app here.

Are there any workarounds to avoid this behaviour? It's easy to accidentally double-click within the Table of Contents, causing the map to zoom unexpectedly.


texts: Why does clipping a line layer using a polygon result in an empty shapefile?

texts:
I am using QGIS 2.2.0 and want to clip a line layer (roads) so that only those inside a polygon are shown. When I do this using Vector->Geoprocessing Tools->Clip I only get an empty shapefile as an output. I have tried transforming the original road shape file from multipart to singlepart before clipping and still have the same problem.

I've included a screenshot of what I am trying to do.

Does anyone have an idea of what is going wrong or something to do differently?


texts: openlayers interaction only on left mouse

texts:
Using openlayers (latest - 4), I have a drawing interaction that draws a circle -

let drawCircle = new ol.interaction.Draw({source: source, type: "Circle"});
map.addInteraction(drawCircle);

Both left mouse and right mouse click have the same behavior with this interaction- they start and stop the drawing of the circle.

I also have a context menu, so I do not want the right click to do anything with the interaction, just the left click.

How do I make this happen?


texts:
Does anyone please knows where to get global surface temperature maps, yearly?
I need separate map for each year from 2002 to 2012 shows the average surface temperature over the globe, would be great if it is in GIS/Envi format.


texts: GML file provides wrong coordinate system

texts:
I have the following field 51A.
It is described by the following GML excerpt:

<gml:featureMember>
    <ogr:av_export_temp.ihkulroo fid="av_export_temp.ihkulroo.183">
      <ogr:geometryProperty><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-63613.5803279425,6594849.11111532 -63338.4982627596,6596691.79998076 -61065.5082310117,6596353.15340253 -58792.3850046138,6596015.89007513 -56519.1291541683,6595680.00990826 -54245.7412499302,6595345.51281189 -51972.2218618083,6595012.39869636 -52241.6978222727,6593169.17366796 -54516.3382854795,6593502.39409221 -56790.8473374676,6593836.99795894 -59065.2244085927,6594172.98535813 -61339.468928864,6594510.35638005 -63613.5803279425,6594849.11111532</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
      <ogr:Copyright>Fishdirectory</ogr:Copyright>
      <ogr:Zone>51A</ogr:Zone>
      <ogr:County>Rogaland</ogr:County>
      <ogr:Timereport1>2018/09/11</ogr: Timereport1 >
      <ogr:Timereport2>2019/08/31</ogr:Timereport2>
    </ogr:av_export_temp.ihkulroo>
  </gml:featureMember>
  <gml:featureMember>

When I looked into GML file and QGIS CSR properties, it says that it uses EPSG:25832 ETRS89 / UTM zone 32N.

If I transform the first coordinate point (-63613.5803279425,6594849.11111532) and then display it on the Google map, it is somewhere in the North Sea, but it should be somewhere around (59.130778, 5.241507) (taken from Google Maps). What am I missing here?


texts: Calculate cloud free median value for landsat 5, 7, 8 image collections in GoogleEarthEngine

texts:
I have an image collection of Landsat 5, Landsat 7 and Landsat 8 images. I want to calculate the median value of my image collection. Before that I want to exclude all cloudy pixels of my collection, so that they are not considered to calculate the median.

The standard function maskL8sr(image) doesn't seem to work. It shows for my median value image clouds inside the image as shown in my images.

First Image: Single Landsat 8 image of the month October

Second Image: Median value of image collection for the month's May till October

My code looks so far like this:

var geoshape = ee.FeatureCollection('ROI');

var col5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR").select(['B1','B2','B3','B4','B5','B7','pixel_qa']);

var col7 = ee.ImageCollection("LANDSAT/LE07/C01/T1_SR").select(['B1','B2','B3','B4','B5','B7','pixel_qa']);

var col8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").select(['B2','B3','B4','B5','B6','B7','pixel_qa'],['B1','B2','B3','B4','B5','B7','pixel_qa']);

var col57 = col5.merge(col7);
var col = col5.merge(col8);

function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var dataset = ee.ImageCollection(col)
          .filterDate('2016-05-01', '2016-10-31')
          .filterBounds(ROI)
          .filter(ee.Filter.eq('WRS_ROW', 52))
          .filter(ee.Filter.lt('CLOUD_COVER', 25))
          .map(maskL8sr);

var medianpixels = dataset.mean();
var medianpixelsclipped = medianpixels.clip(geoshape);


texts: WFS Service CRS error, swapped coordinates?

texts:
I'm having trouble getting the USDA WFS services described here to display features in QGIS.
I've used various versions of this as my WFS service URL: https://sdmdataaccess.nrcs.usda.gov/Spatial/SDMWGS84Geographic.wfs?VERSION=1.1.0 and keep getting this error:

  Layer mapunitpoly: Download of features for layer mapunitpoly failed or partially failed: Download of features failed: Error transferring https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMWGS84GEOGRAPHIC.wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&TYPENAME=mapunitpoly&SRSNAME=urn:ogc:def:crs:EPSG::4326&BBOX=61.66538482856496017,-148.88267110247048208,61.68470628768040598,-148.84056295518394109,urn:ogc:def:crs:EPSG::4326 - server replied: Bad Request. You may attempt reloading the layer with F5

I checked this URL in my browser and get an error about using the correct CRS:

  Error: parameter "srsname" requires a value from the list ("EPSG:4326","EPSG:3857","EPSG:42659")

The error seemed similar to errors prior to the addition of WFS v.1.1, but I am using QGIS 3.10.1-A.

I tried disecting the url and get an error about bounds when I get to this step:

  https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMWGS84GEOGRAPHIC.wfs?SERVICE=WFS&REQUEST=GetFeature&TYPENAME=mapunitpoly&BBOX=61.66538482856496017,-148.88267110247048208,61.68470628768040598,-148.84056295518394109

It appears the syntax has Lat/Long switched?

  Error getting BBOX extent area from WFS request string 24201: Latitude values must be between -90 and 90 degrees.

Using this URL in my browser, I get data back from the server. Yay!

  https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMWGS84GEOGRAPHIC.wfs?SERVICE=WFS&REQUEST=GetFeature&TYPENAME=mapunitpoly&BBOX=-148.88267110247048208,61.66538482856496017,-148.8405629551839410,61.684706287680405989

However, when I try to directly enter this as my URL in the QGIS WFS list, it becomes:

  https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMWGS84GEOGRAPHIC.wfs?BBOX=-148.88267110247048208,61.66538482856496017,-148.8405629551839410,61.684706287680405989&SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&TYPENAME=mapunitpoly&SRSNAME=urn:ogc:def:crs:EPSG::4326&BBOX=60.48852587202364361,-153.37828957730314983,61.6972444843330976,-150.43631069033008885,urn:ogc:def:crs:EPSG::4326

Which returns the same error as at the start. What am I doing wrong?


texts: Convert raster LST 1 km to 250m (R resample)

texts:
I have two files, where the NDVI is at 250m and the LST is at 1k, I want to use resample of R to convert the LST to 250m.

ndvi <- raster("MOD13Q1.A2014017.tif")
lst <- raster("MOD11A2.A2014017.tif")
LST250 <- resample(lst, ndvi, method="ngb", filename="lst_250.tif")

But I get the following error, they could help me

[Edit] Unfortunately, I have discovered that it is not possible to do a function with resample, apparently it can only be called once, otherwise it gives me the error I have shown

resample <- function(){
    ndvi <- raster("MOD13Q1.A2014017.tif")
    lst <- raster("MOD11A2.A2014017.tif")
    LST250 <- resample(lst, ndvi, method="ngb", filename="lst_250.tif")
}

resample()


texts: Why do all rows go missing from File Geodatabase Feature Class to *.dbf conversion?

texts:
I tried to export a File Geodatabase Feature Class to a *.dbf table/shapefile, but there are no records inside. I believe there is no mistake in exporting--it is simple. The doubt is the file geodatabase format--I cannot use it as a shapefile or export it as a shapefile with intact data. Attached is the original data and the exported data.

I am using ArcGIS 10. The data are from the United States census.

I am only interested in one data column/Field from each state. Is there a fast way to open, export, or convert the geodatabase format to *.dbf?


texts: Specify the adminsitrative level in QuickOSM

texts: I am using QuickOSM to download the boundaries of cities into QGIS before saving each city data as a GeoJSON. QuickOSM enables the Key (in my case, 'boundary') and Value ('administrative) to be specified, but I don't know how to specify the administrative level (admin_level=x). However, I need to do this to get the boundaries for the actual city. At present, I get boundaries for the wrong administrative level. For instance, running the query for 'Amsterdam' outputs the boundaries for the Netherlands as a whole. What can I include in the Overpass Query to specify the administrative level?


texts:
I am working with data that I need to be able to display by day gathered, and it's gathered over a period of years.  I would like the user to be able to select a date from a calendar-type picker, and then only display points from that date.

My gut tells me that this must have been done before, but I can't find any examples anywhere.  I'm tempted to dive in and just start coding stuff and making it work, but I'd like to avoid reinventing this particular wheel.  Can anyone point me in a good direction?


texts:
I would like to add a HTML label above each marker in a leaflet map.

Is this possible? Example:

Code for labels:

HTML

CSS

table{
    border-collapse: collapse;
}
td{
    padding:20px;
    border:5px solid black;
}
.red{
    background-color:#F15E66;
}
.yellow{
    background-color:#FFDB64;
}
.orange{
    background-color:#F58326;
}
.blue{
    background-color:#85B1DE;
}


texts:
I'm using OpenLayers 2.13 with GeoExt and mapfish print module for GeoServer.
I've a map composed of a WMS layer and a vector layer. On my web browser everything is fine but when I print this map into a PDF file with mapfish print, layers are shifted and I can't figure out why. Any ideas?

How the map is rendered by OpenLayers/GeoExt

How the map is rendered in PDF file generated by Mapfish print

Part of my code :
Map creation part:

map = new OpenLayers.Map("map", options);
var layer_dpt = new OpenLayers.Layer.WMS("dpt", "linkToGeoserver",
    {
        layers:'base_layer'
    },
    {
        singleTile: 'true'
    }
);
map.addLayer(layer_dpt);
layer_comm = new OpenLayers.Layer.Vector( 'communes',{
strategies: [new OpenLayers.Strategy.Fixed()],
protocol: new OpenLayers.Protocol.HTTP({
url: 'GeoServer WFS URL',
format: new OpenLayers.Format.GeoJSON()
})
});
map.addLayer(layer_comm);

Printing part:

var printProvider = new GeoExt.data.PrintProvider({
    method: "POST", // "POST" recommended for production use
    capabilities: printCapabilities // from the info.json script in the html
});
printProvider.setDpi(printProvider.dpis.getAt(2));
// Our print page. Tells the PrintProvider about the scale and center of
// our page.
printPage = new GeoExt.data.PrintPage({
    printProvider: printProvider,
    customParams: {
        mapTitle: "",
        comment: ""
    }
});
mapPanel = new GeoExt.MapPanel({
    region: "center",
    map: map,
    center: [map.getCenter().lon, map.getCenter().lat],
    zoom: map.zoom
});
new Ext.Panel({
    renderTo: "map",
    layout: "border",
    width: 1200,
    height: 950,
    items: [mapPanel],
    bbar: ["->", {
        text: "Print",
        id: "imprimer",
        handler: function() {
            // convenient way to fit the print page to the visible map area
            var opt = {"mode":"closest"};
            printPage.fit(mapPanel, opt);
            // print the page, optionally including the legend
            printProvider.print(mapPanel, printPage, false);
        }
    }]
});


texts: QGIS Python Script - print "HELLO QGIS!" using a Windows Batch File (.bat)

texts:
I have a python script (as given below) in QGIS which will just print "Hello QGIS!".

import sys
from qgis.core import *

# Initialize QGIS Application
QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
app = QgsApplication([], True)
QgsApplication.initQgis()

# Add the path to Processing framework
sys.path.append('C:/OSGeo4W/apps/qgis/plugins')

# Import and initialize Processing framework
from processing.core.Processing import Processing
Processing.initialize()
import processing

print 'Hello QGIS!'

Now, I want to execute this python script without opening QGIS by using a Windows Batch File script. The .bat file which I created for this is given below.

REM Change OSGEO4W_ROOT to point to the base install folder
set QGIS_PREFIX_PATH=C:\OSGeo4W\apps\qgis
REM Gdal Setup
set GDAL_DATA=C:\OSGeo4W\share\gdal\
REM Python Setup
set PATH=C:\OSGeo4W\bin;C:\OSGeo4W\apps\qgis\bin;%PATH%
SET PYTHONHOME=C:\OSGeo4W\apps\Python27
set PYTHONPATH=C:\OSGeo4W\apps\qgis\python;%PYTHONPATH%

REM Launch python job
python C:/Users/Sreeraj/.qgis2/processing/scripts/hello.py
pause

Unfortunately, I am getting an error :

File "C:/Users/Sreeraj/.qgis2/processing/scripts/hello.py", line 13, in <module>
from processing.core.Processing import Processing
ImportError: No module named processing.core.Processing

I also tried modifying the .bat file as given below.

SET QGIS_PREFIX_PATH=C:\OSGeo4W\apps\qgis
call "C:/OSGeo4W/bin/o4w_env.bat"
SET PYTHONPATH=C:\OSGeo4W\apps\qgis\python;%PYTHONPATH%
SET PATH=%PATH%;C:\OSGeo4W\apps\qgis\bin

cmd /c python "C:/Users/Sreeraj/.qgis2/processing/scripts/hello.py"

Again, I am getting the same error :

ImportError: No module named processing.core.Processing

But, clearly, I have set everything correctly and the processing folder is located in C:/Users/Sreeraj/.qgis2/ .

So, what kind of modification do I have to make in order to run this basic python script correctly using a .bat file ?


texts: Color gradients to show direction in polylines

texts:
I am referring to this very useful hack to change the color gradient of buffered polygons in order to show direction (with a gradient from red to green, for example)

hack

The only problem is that data defined properties is no longer a button in QGIS.  I am running 2.12.  How do I get around this?

Here is an example of the lines that I want to fill with gradient color to show their direction:

Any other suggestions on how to show direction would be welcomed!

This is one of my attempts but actually the color gradients are random


texts: Increasing display performance of ArcGIS Online Basemaps like Bing in ArcMap?

texts:
I'd like to know if there are ways (within ArcGIS itself) to increase the display performance of basemaps within ArcMap. The internet connection was recently upgraded, so it is much faster than before, and by that I mean instead of taking 1 minute to draw the Bing Imagery basemap, it now takes about 30 seconds.

I normally pan and zoom with the basemap off, and once I think I'm in the location I want, I turn the basemap on and take my hands away from the mouse (trying to do anything else guarantees me a spinning circle, and most likely a crash). From this post on the ArcGIS forums, I've changed the location of my display cache from My Documents to my C drive, and I could see an improvement.

I then changed the location to our server, and it now takes about 5 seconds to draw when I pan and zoom. While I am quite happy with this, does the fact that I am using 10.1 make a difference? I changed the cache locations for my colleagues who are using 10.0, and their performance did not improve at all. I don't know what else to do, as I can't do anything about the antivirus/firewall part (that's under IT's control).

It's a bit frustrating, particularly when the basemaps in ESRI Maps For Office and ArcGIS for AutoCAD draw almost instantly, and I can zoom in without worrying my computer will implode. My current settings in ArcMap options are hardware acceleration enabled, and cache location changed.

I've looked at the help on managing your local display cache, but my settings appear to be optimised already, unless I am missing something.


texts: Leaflet Marker Cluster: refreshCluster function issue

texts: I am trying to use the markers.refreshClusters(); function to refresh the marker icons for a new styled icon I am trying to work with to highlight the selected cluster. Unfortunately, I keep getting the markers.refreshClusters is not a function error. Is there a file I am missing for the marker cluster plug-in? I tried updating to latest revision of plug-in but it seems I have the latest.


texts: Random coordinates need to move little right

texts:
I have amcharts map which shows circle on it. I have few records(nearly 200) which do not have locations. These records(circle) I am showing in the map in Sea so that those records clearly will be visible.
I have tried, but circles I need to have little space between circles and those need little right side
This is my code
result.forEach(element => {
if(element.loc=='NULL')
{
  var r = 1000590/111300, // = 100000 meters
  y0 = 27.990741,
  x0 = 153.556045,

  u = Math.random(),
  v = Math.random(),
  w = r * Math.sqrt(u),
  t = 2.3 * Math.PI * v,
  x = w * Math.cos(t),
  y1 = w * Math.sin(t),
  x1 = x / Math.cos(y0),
  newY = y0 + y1,
  newX = x0 + x1
  element.lat=newY
  element.long=newX
}
});

This is my output

Something like this I need but little right side

How can I fix this ?


texts: Adding new field conditioned by analysis that regards other layers

texts:
I'm searching for a method to add a new field in the point vector that has a specific value of a field contained in a polygon vector with QGIS. Here is an example image:
The violet is the point shapefile to which I want to add the number of the polygon (in black) which contains them.

Someone has any idea?


texts: Selecting nearest lines to point layer in QGIS?

texts:
I am a beginner with GIS software and I'm working on QGIS.

I have a problem joining the spatial layers. I have a layer with objects in one city and another layer with the streets. I need to attach each object to its street.

I try to join attributes by location but I have not been able because it does not take into account the objects that are on the edge of the street (the roundabouts and bridges).

So how can I make a join between a point and a layer of layer lines using the closest distance by QGIS.


texts:
I have a point shapefile (Houses.shp) and various lines shapefiles (SemiRadial.shp, MST.shp, Ring.shp, etc) representing respectively houses of a village and possible microgrid network layout connecting those houses together. I would like to create networks using QNetwork plug-in to analyze different characteristics of each network, such as average shortest path between nodes.

As far as I understood, to create a network with QNetwork you need at least one node layer, which would be my house shapefile, and one polyline layer, which would be any of my lines shapefiles. The problem is that during the creation of the network QGIS crashes (tried that several times, with different polyline layers and also with new point/polyline layers created on purpose to see if my existing layers were the problem).

I noticed that even with the crash, the network directory "my_network" was indeed created, but it's not possible to load it in the project using QNetwork --> Load Network because it does not read it. Does anybody have experience using QNetwork and could maybe give me some hints on what I am doing wrong?

I attach images of the various steps of the process in case I am doing some basic silly mistakes.


texts:
I have an RGB orthoimage that has black area instead of NoData, as pictured here:

When I use the SetNull tool, it does remove the black area, however it then turns it into a grayscale orthoimage, which is not what I want.

I know that there is a way around this by using the Clip (Data Management) tool to clip the raster to itself and in there you can set the NoData value to be the value of the black cells (0). However, I am creating a workshop for students and I would like to use a less counterintuitive method.

There is of course the option of running the NoData tool on each band individually and then combining the three bands together, however I need to have a condition where band1, band2, and band3 all equal 0 to set that as NoData. Being able to do this in ModelBuilder would be ideal, but scripting it using ArcPy would also be good.

I am using ArcGIS 10.2 Desktop with Advanced License.


texts: WFS GetFeature with propertyName get response without geometry

texts:
I have a WFS request send to a GeoServer:
http://localhost:8600/geoserver/aa/ows?service=WFS &version=2.0.0&request=GetFeature&typeNames=aa:bb&outputFormat=json&propertyName=attr1,attr2
I cannot get all features with only attr1 and attr2 value. When requesting every feature's geometry the response is null.
When I send a similar request without propertyName parameter, the geometry is back.
I want to have both the geometry and attr1, attr2.


texts:
How can I convert a GOCI satellite image from Radiance to Reflectance?
I have been searching through the internet, but most of them are related to other satellites. I try to calculate the reflectance using the typical approach, yet again, I found that I cannot find proper value of ESUN.

Is there something obvious I have missed or misunderstanded?


texts: Automatisation of Zonal Statistics with python script

texts:
I am using the Zonal Statistics tool of QGIS 2.18.9 to calculate the mean value of a raster within a grid of rectangles.
When I do it "manually", everything works. But as I have to do this for a lot of grids and raster files (more than 100), I want to automate this.
I tried by different ways :

1)Using Graphical Modeler
It failed because you can't specify what you want to get (min, max,mean ...) so it creates all the column with only zeros inside.

2)I tried with a python script.
So I wrote a script and I run it in the Python console of QGIS.
I tried to different ways :

2) a) I read this article : http://spatialgalaxy.net/2012/01/27/qgis-running-scripts-in-the-python-console/
And I succeeded in running his script. I modified it to get what I want :
But python returns me : "NameError: global name 'QgsVectorLayer' is not defined"

So as the man in the article tells about a "simple" way. I tried it too.

2)b) Still using python, this time I am not creating class.
I just write it directly.
As for 2)a) I tried a simple code that I found there :
How to use QGIS Zonal Stats Plugin from Python Console?
It works. So I changed it to obtain what I want :
it gave me nothing, but when I write "zoneStat.calculateStatistics(None)
" in the Python Console it returns 3 and change nothing in my Vector Layer.
Same thing when I want to run it line by line without python script : when I write in the command it returns me 3 also.

So, to put it in a nutshell :
I try to run a python script to automate the QGSZonalStatistics method.
With classes it does not recognize  QgsVectorLayer (so perhaps I forgot to import something)
With a simple dirty script or only by writing in the command, it returns 3 without further information. I read on a forum that it should return 0 if it created the columns.


texts: Cokriging cross validation variances too small?

texts:
I'm looking into kriging the Meuse dataset in R. I first performed ordinary kriging using copper, using krige.cv to calculate residuals. The distribution of prediction variances given by applying krige to my grid and krige.cv was approximately the same. However, when I try doing cokriging, the cross validation prediction variances are 60-75% smaller than before, whilst the prediction variances given by "predict" are roughly the same as with ordinary kriging. Is there a reason to explain this - you can see the cokriging code below.

library(sp)
library(gstat)
data(meuse)
meuselog<-meuse
meuselog["logcadmium"]<-log1p(meuselog$cadmium)
meuselog["logcopper"]<-log(meuselog$copper)
coordinates(meuselog) = ~x+y
proj4string(meuselog)<-CRS("+init=epsg:28992")

data(meuse.grid)
coordinates(meuse.grid) = ~x+y
proj4string(meuse.grid)<-CRS("+init=epsg:28992")
gridded(meuse.grid)=TRUE

lcu.vgm <- variogram(logcopper~1, meuselog)
lcu.fit <- fit.variogram(lcu.vgm, model = vgm(0.5,"Sph",900,0.1))

#Cokriging - copper and cadmium
coppercd<- gstat(NULL, id = "logcopper", form = logcopper ~ 1, data=meuselog)
coppercd<- gstat(coppercd, id = "logcadmium", form = logcadmium ~ 1, data=meuselog)
coppercd.cross<-variogram(coppercd)
plot(coppercd.cross)
# Add variogram models to the gstat object and fit them using LMC
coppercd<-gstat(coppercd, id = "logcopper", model = lcu.fit, fill.all=T)
coppercd<-fit.lmc(coppercd.cross,coppercd)

coppercd.ck <- predict(coppercd, meuse.grid)
set.seed(1234)
coppercd.cv <- gstat.cv(coppercd,nfold=10)
summary(coppercd.ck)
summary(coppercd.cv)


texts: SyntaxError: can't assign to function call, arcpy set parameter as text

texts:
I have a python script loaded into arcgis ModelBuilder, and when I run the model I keep getting an error when the script tries to run, at the line where I am setting the output as a text parameter. The error is SyntaxError: can't assign to function call. I've researched this error but still haven't figured out what the problem is? Here is my script:

import arcpy

mapFrames = r'C:\temp\temp.gdb\D04MapIndex_Master_MapSheetIndex'
fields = ['linkdir','mapdecade']
decadesText = r'P:\Records\GIS\__Databases__\MapIndex\decades_join_table.txt'
arcpy.SetParameterAsText(0) = decadesText

with open (decadesText,'w') as f:
    f.write('dir;decade\n')
    for row in arcpy.da.SearchCursor(mapFrames,fields):
        if row[1] is None or row[0] is None :
            print 'No decade or dir value, skipping'
            arcpy.AddMessage('No decade or dir value, skipping')
        else:
            print 'writing out',row[0],', decade value:',row[1]
            arcpy.AddMessage('writing out '+row[0]+', decade value:'+str(row[1]))
            f.write(row[0]+';'+str(row[1])+'\n')
f.close()


texts: Getting SearchWidget result to JavaScript variable?

texts:
I'm using the following searchWidget and it works fine, I need to know how can I get location details, once I type something and hit enter.
Lets I search for New York and hit enter, then search data should be taken into a variable (search value, location details, etc).

  var searchWidget = new Search({
        sources:[{
          countryCode:"LKR",
          placeholder: "Search",
        }],
        view: this.view,
        includeDefaultSources: false,
        container : document.createElement("div")
      });


texts: Does ArcGIS Collector capture device metadata?

texts: Does ArcGIS Collector capture device metadata along with observed data ie. actual device location, system date, GPS accuracy?


texts:
I'm using OpenLayers 2.13 and Postgresql 9.3 for my development activities.
I want to load region data from database and load into openlayers.

Below result i'm getting from St_astext(geom) in postgresql 9.0 and above

  "MULTIPOLYGON(((77.6620335903023 12.8381894855534,77.6620330917475
  12.8384550503658,77.662583464913 12.8384554346975,77.6625831247336 12.8381894684868,77.66234257162 12.8381895285991,77.6622741102056 12.8381895459662,77.6620335903023 12.8381894855534)))"

But in postgresql 9.3 i'm getting the below result with height parameter while using St_astext(geom)

  "MULTIPOLYGON ZM (((77.6620335903023 12.8381894855534 0 0,77.6620330917475 12.8384550503658 0 0,77.662583464913 12.8384554346975 0 0,77.6625831247336 12.8381894684868 0 0,77.66234257162 12.8381895285991 0 0,77.6622741102056 12.8381895459662 0 0,77.662033590 (...)"

Below is the OpenLayers code to featch the wkt format

  var feature = new OpenLayers.Feature.Vector(new
  OpenLayers.Geometry.fromWKT("MULTIPOLYGON ZM (((77.6620330345526
  12.8382017120463 0 0,77.6620403164933 12.8384556061155 0 0,77.6625840206627 12.8384543231982 0 0,77.6625747884884
  12.8381889127371 0 0,77.6620330345526 12.8382017120463 0 0)))"), null,{

  label : "TEST",
  irregular: true,
  cursor:'pointer',
  fontSize: "14px",
  fontFamily: "arial, monospace",
  fontWeight: "bolder",
  labelAlign: "cm",
  fillColor: "#9fdaef",
  fillOpacity: 0.3,
  strokeWidth:4,
  strokeColor: "#005596"
  });

I'm getting an error like below

  TypeError: a.geometry.getBounds is not a function in OpenLayers.js
  (line 611, col 334)

By using Multipolygon without height parameter i can able to load region data in Openlayers

but, with the height parameter i can't.

Please guide me to load wkt value with height parameter in Openlayers.


texts: Can a custom add-in execute code upon install of the .esriaddin

texts:
I have made a custom add-in using C#.NET(a toolbar that contains buttons) and I was wondering if it is possible to have some code execute upon the install of the .esriaddin file. Basically I want to create a directory with some map templates on the users computer. I can just add the

DirectoryInfo di = Directory.CreateDirectory(@"Path"); etc....

to the on click event of the buttons and have it only execute if the folder doesn't exist but I was wondering if there was a way to execute it on install so I don't have to include it within the buttons.

Edit: I have looked at this and it is similar though I'm not using a custom installer.


texts: SQLite create view syntax error in ArcPy Python script tool only

texts:
I'm working with a SQLite database in a python script tool in ArcGIS Desktop 10.4. When I try to create a view on the database in SQLite Studio (using the sql_file pasted below) or from the command line (using the code and the sql_file below), the SQL statement works fine. When running from a Python script tool (using the code and the sql_file below), I get a syntax error ("near WITH: syntax error").

Here is the code that executes the sql statement:

import sqlite3
from sqlite3 import Error
import sys

with open(sql_file, 'r') as file:
        sql = file.read()
    try:
        c = conn.cursor()
        c.execute(sql)
    except Error as e:
        conn.close()
        AddMessage(e)
        sys.exit(0)

The contents of the sql_file are:

CREATE VIEW IF NOT EXISTS view_desktop_results AS
WITH precip_unstack AS (
        SELECT map_unit_id,
               (CASE WHEN precip = 'Arid' THEN proportion ELSE 0 END) AS arid,
               (CASE WHEN precip = 'Mesic' THEN proportion ELSE 0 END) AS mesic
          FROM project_precip
    ),
    mgmt_cat_unstack AS (
        SELECT map_unit_id,
               (CASE WHEN mgmt_cat = 'PHMA' THEN proportion ELSE 0 END) AS phma,
               (CASE WHEN mgmt_cat = 'GHMA' THEN proportion ELSE 0 END) AS ghma,
               (CASE WHEN mgmt_cat = 'OHMA' THEN proportion ELSE 0 END) AS ohma
          FROM project_mgmt_cats
    ),
    wmz_unstack AS (
        SELECT map_unit_id,
               (CASE WHEN wmz = 'MZ III' THEN proportion ELSE 0 END) AS mz3,
               (CASE WHEN wmz = 'MZ IV' THEN proportion ELSE 0 END) AS mz4,
               (CASE WHEN wmz = 'MZ V' THEN proportion ELSE 0 END) AS mz5
          FROM project_wmz
    )
    SELECT m.map_unit_id,
           m.map_unit_name,
           m.meadow,
           (CASE WHEN m.conifer_phase IS NULL THEN 'None' END) AS conifer_phase,
           m.indirect_benefits_area,
           m.map_unit_area,
           m.no_transects,
           m.spring_hsi,
           m.summer_hsi,
           m.winter_hsi,
           cl.ls_breed as current_breed,
           cl.ls_summer as current_summer,
           cl.ls_winter as current_winter,
           pl.ls_breed as projected_breed,
           pl.ls_summer as projected_summer,
           pl.ls_winter as projected_winter,
           mgmt.phma,
           mgmt.ghma,
           mgmt.ohma,
           wmz.mz3,
           wmz.mz4,
           wmz.mz5,
           precip.arid,
           precip.mesic
      FROM map_units AS m
           LEFT JOIN
           current_ls AS cl ON m.map_unit_id = cl.map_unit_id
           LEFT JOIN
           projected_ls AS pl ON m.map_unit_id = pl.map_unit_id
           LEFT JOIN
           mgmt_cat_unstack AS mgmt ON m.map_unit_id = mgmt.map_unit_id
           LEFT JOIN
           wmz_unstack AS wmz ON m.map_unit_id = wmz.map_unit_id
           LEFT JOIN
           precip_unstack AS precip ON m.map_unit_id = precip.map_unit_id;

Note that I'm running a bunch of other SQL statements to create the tables referenced in the SQL statement above. Creating tables and inserting rows all works fine. Only views seem to be affected. Again, this statement works outside of ArcGIS.

I'm thinking this may be a version issue between what is distributed with ArcGIS Desktop and the more recent SQLite distributions I'm using outside of the ArcGIS environment.

On ArcGIS Desktop 10.4 Basic, running Python version: 2.7.10 (default, May 23 2015, 09:40:32) [MSC v.1500 32 bit (Intel)] with Windows 8. sqlite3.version has given me 2.6.0 and sqlite3.sqlite_version has given 3.6.21


texts: How to consider more than one no data value in utility gdalwarp of GDAL

texts: I have to do a bilinear interpolation on an image with two no data values, using the utility gdalwarp of GDAL (http://www.gdal.org/gdalwarp.html), so these two values shouldn't be considered in interpolation. How can I do this? I evaluated the gdalwarp options "-srcnodata" and "-dstnodata", but they specify several no data values, with each value corresponding to one band. The case of my image is different, because I have only one band with multiple no data values.


texts: SQL Server ArcSDE ST_Raster Installation and Configuration Error ".dll"

texts:
I have successfully installed ST_Raster Library in SQL Server machine following this tutorial

http://desktop.arcgis.com/en/arcmap/10.4/manage-data/using-sql-with-gdbs/install-st-raster-oracle.htm

I have used ArcGIS 10.4 Configuration wizard for import and export of configurations
After installation i have also run this query.

  SELECT sde.st_raster_util_getVersion();

After installation i have also run this query.

  SELECT sde.st_raster_util_getVersion();

This query successfully returned the "1000" that means i have successfully installed the ST_raster, i am also able to create a table in SQL Server and from design of the table i can see the type is set ot ST_Raster.

  CREATE TABLE sde.urban_area_sde (raster sde.ST_Raster, name
  varchar(255));

When initializing the document as tutorial say i used below query to initialize the raster column

EXEC sde.ST_Raster_Util_Initialize 'AgriAhDB','sde','urban_area_sde','raster',4326,NULL,'DEFAULTS'

It was returned as error.

  Msg 6522, Level 16, State 1, Procedure st_raster_util_initialize, Line
  0 [Batch Start Line 6] A .NET Framework error occurred during
  execution of user-defined routine or aggregate
  "st_raster_util_initialize":  System.DllNotFoundException: Unable to
  load DLL 'pe.dll': The specified module could not be found. (Exception
  from HRESULT: 0x8007007E) System.DllNotFoundException:     at
  ESRI_ArcSDE.STRasterUtil.SdeRasterAPI.pe_factory_coordsys(Int32
  projectionID)    at
  ESRI_ArcSDE.STRasterUtil.StoredProcedures.getCoordSysFromPID(SqlInt32
  projectionID, String& srtext)    at
  ESRI_ArcSDE.STRasterUtil.StoredProcedures.getCoordSys(SqlInt32
  projectionID, SqlString projectionFile, String& srtext)    at
  ESRI_ArcSDE.STRasterUtil.StoredProcedures.createRasterColumn(SqlString
  database, SqlString owner, SqlString tableName, SqlString columnName,
  SqlInt32 projectionID, SqlString projectionFile, SqlString keyword,
  SqlInt32& rastercolumnID, SqlInt32& regID)    at
  ESRI_ArcSDE.STRasterUtil.StoredProcedures.st_raster_util_initialize(SqlString
  database, SqlString owner, SqlString tableName, SqlString columnName,
  SqlInt32 projectionID, SqlString projectionFile, SqlString keyword) .

I have restarted the SQL Server service, but it was unsuccessful in resolving the error. How can I resolve this error,? Is there any specified library that I need to place to use the ST_Raster procedure?


texts:
I am working with some rasters that have been produced with a wrong projection.
The Spatial Reference is saved as WGS84 (EPSG:4326) but the unit are clearly in meters and not degrees. I have tried both with ArcMap and QGIS.

If I try to project to UTM, something happens but the units are clearly degrees.

Is there a way to force, or manually overwrite the projection?

Unfortunately there's no chance to get the correct data from the source that sent them to me.


texts: ArcGIS JavaScript setRequestPreCallback is partially working?

texts:
I have an api sitting in front of AGS that requires it's own authentication. I am using setRequestPreCallback as stated by the ESRI docs.

In the image below, you can see that the first 2 requests are returning 200. The first request is CORS Preflight and the second request is the JS api getting the json info from the MapServer.

The second request uses the function I have set in setRequestPreCallback(...) to properly set a header value with the token. However, whenever an image is requested from the MapServer (via export), the function I have in setRequestPreCallback is not called and the token is not set, leaving me with a bunch of 401 Unauthorized requests.

I am running AGS 10.3 with JS API 3.14

Here is the JS I am using:

require([
    "esri/map",
    'esri/request',
    'esri/config',
    "esri/layers/ArcGISDynamicMapServiceLayer",
    "dojo/domReady!"
], function(Map, esriRequest, esriConfig, ArcGISDynamicMapServiceLayer) {

  esriConfig.defaults.io.corsEnabledServers.push("my.url.com");

  var map = new Map("mapDiv", {
    center: [-118, 34.5],
    zoom: 8,
    basemap: "topo"
  });

  function myCallbackFunction(args) {

    if (args.url.indexOf('my.url.com') === -1) {
        return args;
    }

    args.headers = args.headers || {};
    args.headers.token = 'token';

    return args;
  }

  esriRequest.setRequestPreCallback(myCallbackFunction);

  map.addLayer(new ArcGISDynamicMapServiceLayer("https://my.url.com/map/Assets_Service/MapServer"));
});

Any ideas?


texts: Using regex to find an attribute with no letters in FME Workbench

texts:
I need to divide address elements to lines. I am using 'TestFilter' transformator in FME Workbench to separate attributes with no letters. For this purpose I am using IF statement and regular expressions.

This is how Test Clauses look like:

Left Value        Operator           Right Value
HOUSE_NAME        Matches Regex      @FindRegEx(@Value(HOUSE_NUMBER),^[^A-Z]+$)

I am not sure if the syntax is right.


texts:
I have a set of points at which I want to find the minimum spanning tree with PostGIS. I have no lines between them, I only have the starting point from which to start building the tree (the Multilinestring)
I don't know how to start coding this, would it be better to do it with a recursive query? maybe implement the Prim or Kruskal algorithms using the distances between the points?
For now, I have a table with the points (id, geom) and the starting point (start_point: = getStartPoint (points))

MST: minimum weight spanning tree is a subset of the edges of a connected, edge-weighted undirected graph that connects all the vertices together, without any cycles and with the minimum possible total edge weight.


texts: How or where to install Python plugins in OSX? (QGIS 2.18)

texts:
For a school assignment we are using QGIS 2.18 (professor's suggestion) and 2 plugins he gave us. I am working on an iMac. The short manual he gave us is for use of QGIS under Windows.

I have QGIS 2.18 installed, now I need to put the plugins in a place the program knows where to look for. Does anybody have a clue how I can do this?


texts: Circle created using ST_Buffer displays elongated on map

texts:
I am buffering a point to create a circle in PostGIS like:

ST_Buffer(geo::geometry, .0002, 8)

However, when displaying the polygon in Bing maps it ends up a vertically oriented ellipse or elongated circle.

Since both the Geography datatype and Bing both use the WGS84 datum, this is unexpected.


texts:
Could you please help me to convert OpenLayers (srid 3943) extend to Google bounds. I want to geocode from address in my map area so I do this:

var googleProjection = new OpenLayers.Projection("EPSG:900913");
var geocoder = new google.maps.Geocoder();
var address = document.getElementById("geocodeAddress").value;
popupGeo.destroy();

var northEast = new OpenLayers.LonLat(map.maxExtent.right,map.maxExtent.top);
var southWest = new OpenLayers.LonLat(map.maxExtent.left, map.maxExtent.bottom);
var convertedNorthEast = new northEast.transform(map.getProjectionObject(), googleProjection);
var convertedSouthWest = new southWest.transform(map.getProjectionObject(), googleProjection);
var googleBounds = new google.maps.LatLngBounds(convertedSouthWest, convertedNorthEast);

//var convertedBounds = map.maxExtent.transform(map.getProjectionObject(), googleProjection);
//var googleBounds = new google.maps.LatLngBounds(
//    new google.maps.LatLng(convertedBounds.bottom, convertedBounds.left),
//    new google.maps.LatLng(convertedBounds.top, convertedBounds.right));//, 'bounds': googleBounds

geocoder.geocode({ 'address': address }, function (results, status) {

    if (status == google.maps.GeocoderStatus.OK) {

        var newCenter4326 = new OpenLayers.LonLat(results[0].geometry.location.lng(), results[0].geometry.location.lat());
        var newCenter = newCenter4326.transform(googleProjection, map.getProjectionObject());

        map.setCenter(newCenter, GPSLocationZoomLevel);

    } else {
        alert("Impossible de déterminer la position : " + status);
    }
});

but I get an error Phi2z no convergence when I try to convert OpenLayers corners to Google epsg4326


texts:
I have an OpenLayers map with point clustering, but I want to apply a much finer granularity to the data represented in the map.  Specifically, I want to stratify the clusters according to an attribute on the feature (we'll call it 'category').

So, say I have five categories: effectively I need to cluster points only with other points from the same category.  I suppose this could be accomplished with five different layers, but this quickly becomes cumbersome as the number of categories gets large (I have 50+ at the moment).

Has anyone seen anything like this with OpenLayers?


texts: Is there a QGIS plugin that allows users to interface with a PHP form page?

texts:
I've been looking for several days for a QGIS plugin that allows users to interface with a PHP form page.

Does anyone know something that works in this way?


texts: GeoServer ImageMosaic PostGIS index update

texts:
I use GeoServer where I created an ImageMosaic storage following this manual.
When I create a new datastore, all *.tif files are indexed into the DB (PostgreSQL with PostGIS). How can I create new index entry in the DB, when I put new geotiff into my data folder, or reindex all geotiff files in the DB?

My config files:

datastore.properties:

SPI=org.geotools.data.postgis.PostgisNGDataStoreFactory
host=localhost
port=5432
database=GeoServerRaster
schema=public
user=geoserver
passwd=*******
Loose\ bbox=true
Estimated\ extends=false
validate\ connections=true
Connection\ timeout=10
preparedStatements=true

indexer.properties:

Caching=false
TimeAttribute=ingestion
ElevationAttribute=elevation
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date,elevation:Integer
PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](ingestion),DoubleFileNameExtractorSPI[elevationregex](elevation)

timeregex.properties:

regex=[0-9]{8}T[0-9]{6}Z(\?!.\*[0-9]{8}T[0-9]{6}Z.\*)

elevationregex.properties:

regex=(?<=_)(\\d{1,5})(?=_)

files in directory:

tc_10000_20140807T033115Z.tif
tc_10000_20140807T040046Z.tif
tc_10000_20140807T043115Z.tif
...


texts: Unwanted spaces being placed in shapefile text fields by ArcGIS Desktop?

texts:
I have a table with a row Page_Number_N some of the rows have values and others just have a space. For example record 4 is just a space. Field values are Strings.

1
2
3

5

When I start an edit session and delete the space it deletes it. However, as soon as I save edits, the space in the data returns.

I have tried removing the space in the field calc without edit session turned on and the result is the same.

I need this record to equal "" ie. No spaces or null and I can't seem to get ArcGIS Desktop to cooperate.

Does anyone know how to fix this?


texts: Geolocation Ptojection Transform is not working for custom coordinate system

texts:
I am using Openlayers API and I have a custom projection and I defined it to my geoserver named "EPSG:51100".

I have a custom tile map.

        var tms2013 = new OpenLayers.Layer.TMS("2013", "http://mydomin/uydu/2013/", {
                type: 'png',
                getURL: get_my_url,
                isBaseLayer: true,
                tileOrigin: tileOrigin,
                serverResolutions: [.......]
            }
        );

I have generated a geolocation service to zoom on map.

        var geolocation = new OpenLayers.Control.Geolocate({
            bind:false,
            geolocationOptions:{
                enableHighAccuracy: false,
                maximumAge: 0,
                timeout: 7000
            }
        });
        geolocation.events.register("locationupdated", geolocation, onLocationUpdated);

        function onLocationUpdated(e){
            var inProjection =  new OpenLayers.Projection("EPSG:4326");
            var outProjection = map.getProjectionObject();
            var position = new OpenLayers.LonLat(
                 e.position.coords.longitude, e.position.coords.latitude);
            position.transform(inProjection, outProjection);

            map.setCenter(position, 9);
        }

I want to transform the incoming coordinate from geolocate object to my coordinate system but not working.

incoming coordinages like this:
lat: 60.763678
lon: 59.896539799999996

after transofrom process :
lat: 60.763678
lon: 59.8965398

But it should be like this.
lat: 60763.6780
lon: 5989653.9870

I mean the transformation is not working.


texts:
I have these two polygons:

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

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                    55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                    55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                                   242), variable2 = c(235, 464), row.names = c("p1", "p2")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf1, col="red")

coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                     55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                      variable2 = c(3), row.names = c("p1a")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf2, col="yellow", add=T)

I want to cut out parts of spdf1 that are intersected by spdf2. However, I want spdf1 to remain as a SpatialPolygonsDataFrame and to retain any information contained within spdf1@data.

I¹ve tried gDifference as follows, which cuts out parts of spdf1 that are
intersected by spdf2, but then converts spdf1 to SpatialPolygons (i.e. discarding information contained in spdf1@data).

gDifference(spdf1, spdf2, byid=T)

How can I cut into spdf1 with spdf2 but retain data contained in spdf1@data? I have checked and tried this similar question without how to overlay a polygon over SpatialPointsDataFrame and preserving the SPDF data?


texts:
QGIS (1.7) systematically applies the wrong CRS when adding layers (or shapefiles from a given lot). They should be NAD83__MTM9 (EPSG:32189), as can be inferred from the prj-file, but the properties/general/options field show NAD27(CGQ77)/SCPQ_zone9. (EPSG:2015)

Granted from there, I can specify/apply the correct CRS and the layers will line up. But it is a drag to do, and redo again, for many shapefiles. Is there anyway to correct that behavior or a better approach?

Update: indeed corrected with the installation of 1.8


texts:
I face a crucial problem with ODC. I try to do a spatial join with a shapefile (polygon drill). I have followed the following steps:

Read the shapefile and get the first polygon to test
Ensure that projection of the polygon is the same with the polygon of the raw data indexed in the ODC
Make a query based on the geopolygon attribute.

The query does return an empty xarray. However, when I remove the geopolygon attribute from the query it does return all the data. Any idea?
import datacube
import math
# import folium
import calendar
import ipywidgets
import numpy as np
# import geopandas as gpd
import matplotlib as mpl
import matplotlib.cm as cm
from matplotlib import colors as mcolours
import matplotlib.patheffects as PathEffects
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from datetime import datetime
from pyproj import Proj, transform
from datacube.utils import geometry
from datacube.utils.geometry import CRS,Geometry
import shapely
from shapely.geometry import shape
import fiona
import rasterio.features
from fiona.crs import from_epsg
dc = datacube.Datacube(app="test",config="/home/noa/datacube.conf")
shape_file = '/home/noa/Desktop/cyprus36SWDparcels3857.shp'
ds = fiona.open(shape_file)
crs = geometry.CRS(ds.crs_wkt)
for f in ds:
    feature_geom = f['geometry']
    print(feature_geom['extent'])
    geom = geometry.Geometry(geom=feature_geom,crs=crs)
    bounds = shape(feature_geom).bounds
    if 'MULTIPOLYGON' in geom.wkt:
        continue
    break
product = 'sentinel1' # 's2a_sen2cor_granule'
query = {
        'geopolygon': geom,
        'time': ('2018-01-01', '2018-01-31')
}
data = dc.load(product=product, output_crs="EPSG:3857",measurements=['vv'],resolution=(-10,10),**query)


texts: Divide irregularly shaped polygon into 3 equal parts

texts:
I have a polygon object of a transect line that I would like to divide into 3 equal parts that are one-after-the-other along the transect line. So far I have only been able to slice the polygon vertically, which results in non-continuous parts.

Here is a reproducible example using Harvard Forest data, available through the following link

After you have downloaded and unzipped the data, you will see a folder called HARV. In this folder is the shapefile, HARV_roads.shp

Import HARV_roads.shp

lines_HARV <- readOGR("/Your file path/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")

Get footpath lines only

footpath_HARV <- lines_HARV[lines_HARV$TYPE == "footpath",]

Plot the footpath

plot(footpath_HARV,
     lwd=6,
     main="NEON Harvard Forest Field Site\n Footpath")

Create a thicker polygon object from the footpath line using buffer, to simulate a transect line with X width

require(raster)
footpath_buffer <- buffer(footpath_HARV,width=40)

Plot footpath_buffer

plot(footpath_buffer)

Chop the buffer into 3 equal parts using Barry Rowlingson's code (below) from this post

makeVchopper <- function(pol){
  bb = bbox(pol)
  delta = (bb[2,2] - bb[2,1])/10
  xmin = bb[1,1]-delta
  ymin = bb[2,1]-delta
  ymax = bb[2,2]+delta

  choppoly = function(xmax){
    readWKT(sprintf("POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))",
                    xmin,ymin, xmin,ymax, xmax,ymax, xmax,ymin,
                    xmin,ymin))
  }
  choppoly
}

slicer <- function(pol, xmin, xmax){
  bb = bbox(pol)
  delta = (bb[2,2] - bb[2,1])/10
  ymax = bb[2,2] + delta
  ymin = bb[2,1] - delta
  r = readWKT(sprintf("POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))",
                      xmin,ymin, xmin,ymax, xmax,ymax, xmax,ymin, xmin,ymin))
  gIntersection(pol,r)
}

chop_thirds <- function(pol, fractions=c(1/3, 2/3)){
  chopper = makeVchopper(pol)
  bb = bbox(pol)
  xmin = bb[1,1]
  xmax = bb[1,2]

  totalArea = gArea(pol)

  chopped_area = function(x){
    gArea(gIntersection(chopper(x),pol))
  }

  edges = lapply(fractions, function(fraction){
    target = totalArea * fraction
    target_function = function(x){
      chopped_area(x) - target
    }
    uniroot(target_function, lower=xmin, upper=xmax)$root
  })

  xdelta = (xmax-xmin)/10
  chops = matrix(c(xmin-xdelta, rep(edges,rep(2,length(edges))),
                   xmax+xdelta), ncol=2, byrow=TRUE)
  apply(chops, 1, function(edges){
    slicer(pol, edges[1], edges[2])
  })

}

Divide footpath_buffer into 3 equal parts:

parts.footpath <- chop_thirds(footpath_buffer)

Plot footpath_buffer and the 3 equal parts

plot(footpath_buffer)
plot(parts.footpath[[1]], add=TRUE, col=1,border=F)
plot(parts.footpath[[2]], add=TRUE, col=2,border=F)
plot(parts.footpath[[3]], add=TRUE, col=3,border=F)

footpath_buffer is now divided into 3 equal parts, but the parts are broken up. I am hoping to achieve something like this, where each part is continuous along the transect:


texts: Why would programmatically setting the extent of a map on Page_PreRender be really slow?

texts:
We are using the ArcGIS 10 web adf with C# in a web application. We just recently upgraded everything from 9.3 to 10 and moved all the hosted stuff over to a new (more powerful) server. Still on the same network.

During the Page_PreRender event, we get the full extent envelope of our map, which is fast, and then set our map's extent to it.

Code in question...

Envelope fullExtent = AgsHelper.GetExtent(Map1, ExtentType.Map);
Map1.Extent = fullExtent;

The second line there takes a really long time (like 30-40 seconds).

I can't claim to be an ArcGIS expert. I'm wondering if anyone has ideas of what I could look at or what it might be.

I decompiled the code that does the setting of the extent and took a look at the code, but I can't really tell much without being able to debug into it and I can't do that since I don't have a reflector license and my trial is up.

Pretty sure this is some sort of SDE or ArcGIS Server issue, I just don't know what it would be.

Thanks!


texts: OpenLayers3 - Add a source.vector from String, not url (for Layers)

texts:
The current way to add layer from KML file is:

    var urlToKml = someUrl;
    var vector = new ol.layer.Vector({
               source: new ol.source.Vector({
                        url: urlToKml,
                        format: new ol.format.KML()
                    }),
                    visible: true
                });

   map.addLayer(vector);

I'd like to load my source from a string

(ie: where urlToKml = "<kml xmlns='http://www.opengis.net/kml/2.2' xmlns:gx='http://www.google.com/kml/ext/2.2' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xsi:schemaLocation='http://www.opengis.net/kml/2.2 https://developers.google.com/kml/schema/kml22gx.xsd'><Placemark><description>hrt</description><Style><IconStyle><Icon><href>data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAUAAAAFCAYAAACNbyblAAAAHElEQVQI12P4//8/w38GIAXDIBKE0DHxgljNBAAO9TXL0Y4OHwAAAABJRU5ErkJggg==</href><gx:w>5</gx:w><gx:h>5</gx:h></Icon><hotSpot x='2.5' y='2.5' xunits='pixels' yunits='pixels'/></IconStyle></Style><Point><coordinates>2.335711304736855,49.10454920376088</coordinates></Point></Placemark></kml>"; )

Is it possible and how?


texts:
I'm trying to count how many star points there are in 500 meters radius around each property in the month of its transaction. I created buffers around each property transaction firstly.

Since each transaction happened in different months, and there are different distribution of star points in different months, I use Select By Attributes Tool and for example, input "month"=1 for property buffer and "month"=1 for star points respectively, and then I use spatial join tool to get join counts. (Please tell me if this process is wrong)

However, I need to repeat this process for about 10 years, so I'm wondering to know is it possible to let ModelBuilder do it automatically?

I have tried to use Batch for twice, but unfortunately, my computer showed not responding.

I am trying to learn to use  iterators, and i found that every model can only have one iterator. Since I need to change two expressions of select for two input features, I'm wondering to know is it possible to do it using only one iterator?

(I suppose I am still not quiet familiar with how iterators work, but I'm trying to learning it.)


texts: Calculating and using area of polygons in PyQGIS?

texts:
I am using the python code found in Programmatically finding polygons which are >90% overlapped by another vector polygon layer using QGIS?.

I get the following message

error: NameError: name 'area1' is not defined

Here is the concerned part of the code:

for i, feat1 in enumerate(feats_lyr1):
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area1 = feat1.geometry().intersection(feat2.geometry()).area()
            area2 = feat1.geometry().area()
    crit =area1/area2
    if crit > 0.9:
        selected_feats.append(feat1)

After checking lots of things (geometry, CRS, ..), the only thing I see is that i have null area for couple of polygons.

Is that possible that small areas, close to zero, are considered to be zero if considered as integer number?
how to ensure the use of decimal number?


texts: GeoServer CQL cross layer filtering: querySingle/queryCollection for non-geometry attribute return

texts:
I have only been able to find examples where geometry attributes are returned from cross layer filtering functions, but the GeoServer documentation says:

  A common use case is to extract a geometry-valued attribute, and then use the value(s) in a spatial predicate against a geometry attribute in the primary layer.

Which leads me to believe that it was intended to be used to return other data.

I'm looking for examples of how that might be done in a CQL filter.

I have attempted something like this:
http://localhost:8181/geoserver/test/wms?service=WMS&version=1.1.0&request=GetMap&layers=test:layer1,test:layer2&styles=&bbox=-81.54053,41.61621,-77.33826,43.41248&width=768&height=330&srs=EPSG:4326&format=application/openlayers&cql_filter=INCLUDE;property2a=querySingle('test:layer1','property1b','property1a=123')

Where only records in layer2 that have layer2.property2a = layer1.property1b are returned.

I get the exception: "Could not parse CQL filter list."

Similarly:
http://localhost:8181/geoserver/test/wms?service=WMS&version=1.1.0&request=GetMap&layers=test:layer1,test:layer2&styles=&bbox=-81.54053,41.61621,-77.33826,43.41248&width=768&height=330&srs=EPSG:4326&format=application/openlayers&cql_filter=INCLUDE;property2a IN queryCollection('test:layer1','property1b','INCLUDE')

Where the values of layer2 will only be returned if layer2.property2a has a value in the returned values from layer1.property1b


texts: Restrict Feature Info fetched with the WMS GetFeatureInfo request with GeoServer

texts:
I followed a GeoEXT example to create a popup menu that gets WMS GetFeatureInfo on click on the map. What I want to do is to get some specific fields not all the fields; I tried the propertyName (with name1,name2 and name3) and defined the wanted fields but it doesn't work.  I also want to change the title by bringing a specific value not the fid, but it doesn't work, here's the code:

 feature Control, a "button" control
    action = new GeoExt.Action({
    control: new OpenLayers.Control.WMSGetFeatureInfo({
    autoActivate: true,
    infoFormat: "application/vnd.ogc.gml",
    maxFeatures: 6,
    layers: [layers1],
    propertyName: "name1,name2,name3",
    eventListeners: {
    "getfeatureinfo": function(e) {
        var items = [];
        Ext.each(e.features, function(feature) {
            items.push({
               xtype: "propertygrid",
               title: feature.fid, // here I changed it to feature.name1 it doesn't work
               source: feature.attributes
            });
        });
        new GeoExt.Popup({
            title: "infos",
            width: 250,
            height: 300,
            layout: "accordion",
            map: map,
            location: e.xy,
            items: items
        }).show();
      }
    }
   }),
    map: map,
    text: "infos",
    tooltip: "Infos "
      });
    actions["infos"] = action;
    toolbarItems.push(action);
    toolbarItems.push("-");


texts: Exporting table without .geo and system:index in Google Earth Engine

texts:
How do I export a table (CSV) without .geo and system:index in Google Earth Engine?
These two parameters are automatically added when I make a feature.
What I'm trying to do is extract some properties from Image Collection and export CSV of them.
var featColl = someImgCollection.map(function (image) {
    var someValA = image.get("someValueA");
    var someValB = image.get("someValueB");
    var object = {
        someValueA: someValA,
        someValueB: someValB
    }
    var feature = ee.Feature(null, object);
    return feature;
});

Export.table.toDrive({
    collection: featureCollection,
    fileFormat: 'CSV'
});

Then I'd like to remove "system:index" and ".geo" from the result.

system:index
someValueA
someValueB
.geo

2020_01_01
20
34

2020_01_02
14
21


texts:
In QGIS Versions before 3.4 there was a tool called "Extract-Z" values (Not drape!) in the Vector Geometry Toolbox...

Anybody know where it has gone?

I have not found it, neither in the  processing toolbox nor in the vector tools...


texts: Fiona - Preffered method for defining a schema

texts:
I'm fairly new to the world of python and GIS - so this is definitely a beginner question.

I am following Tom Macwrights small guide to Shapely and Fiona, and I noticed that even when opening the previously created shapefile, it doesn't contain the schema I defined using:

schema = {
      'geometry': 'Point',
      'properties': { 'name': 'str' } }

This is evident because when I try to call it using schema = input.schema.copy(), it spits out AttributeError: 'function' object has no attribute 'schema'

Reading the Fiona documentation, it states that:

  the schema of its record type (a vector file has a single type of record, remember) is accessed via a read-only schema attribute.

Does that mean that it is not possible to store a schema in the shapefile,a nd therefore I'm always required set up the variables in the beginning, rather than reading from the file using input.schema.copy()? Or do I setup the schema using Records?

If so, how do you read from a csv file (as in Tom's example) and place the data into the records to be read from input.schema.copy()? is it helpful to do this, or am I just complicating things and should define the schema every time i open the shapefile?

Thanks for the help - I'm really interested in developing my understanding of the processes occurring.


texts: How to delete features/records in ArcGIS Server using Python

texts:
How can I use python to delete records from a table, or features from a featureclass, when they are stored in ArcGIS Server?

using REST "delete features" in arcgis api for javascript shows how to do this in JavaScript - what is the Python syntax?


texts: Getting a flat surface from DSM for height calculation?

texts:
I have a DSM done trough photogrammetry with the software Agisoft Photoscan of a forestry experiment and I need to accomplish an analysis of the mean height from it. My problem is that the profile has a slope and I need a flat surface for the calculations.

This is the exact profile:

I have tried with the Trend tool of ArcGIS but it doesn't work well, still appearing differences. The calculation needs to be very precise because the heights are between 0 to 2.7 meters.

This is one of the results with the Trend tool, and the slope is still appearing. The surface is higher in the top than in the bottom and that doesn't make sense.

Anyone have experience with this?


texts: DBMS Table Not Found with SQLite and ArcGIS 10.3

texts:
I'm trying to set up a spatially enabled SQLite database with ArcGIS 10.3. I've gone through the following steps:

Create a database using the sqlite3 utility: sqlite3 mydatabase.sqlite
Load the relevant extensions: SELECT load_extension('c:\Program Files (x86)\ArcGIS\Desktop10.3\DatabaseSupport\SQLite\Windows32\stgeometry_sqlite.dll','SDE_SQL_funcs_init');
Create the tables to add the ST_Geometry tables to the database: SELECT CreateOGCTables();

After following the above, I can browse to the database in ArcCatalog. I can copy feature classes from personal geodatabases into it using Feature Class to Feature Class. However, whenever I attempt to buffer this feature class in the SQLite database, with the same database the output, I get the error "DBMS table not found".

If a run the same buffer operation from a personal geodatabase to the SQLite database, it succeeds. If I run the same operation from SQLite to a personal geodatabase, it fails.

The database name and feature class names are all lowercase, per the advice of this SE post.

It seems that something is wrong with the data in SQLite, but I have no idea what.


texts: Extracting features from raster data using ArcGIS Desktop?

texts:
I want to extract flood related pixels from my imagery.

I calculated mLSWI(Modified Land Surface Water Index).

According to the literatures, mLSWI >= 0.5 is flood pixels.

How do I identify the flood pixels in my imagery?

I am using ArcGIS Desktop 10.0.


texts:
How can sounding data be used for bathymetry when depth is positive? I have downloaded DEM (digital elevation model) as shown in the image below but need bathymetry data for H. All I have found for H is hydrograph sounding data which can be downloaded in xyz format but data does not match when trying to merge both DEM and sounding data in ArcGIS using spatial analyst/sample tool.
(green elevation less than 0. orange elevation=0)

It is demonstrated how DEM and sounding data does not match.

H05455 XYZ Data


texts: How to add links to a GIS road network, for example, for ferry routes?

texts:
I am using PostgreSQL with PostGIS (for spatial queries like ST_LineLocatePoint()) and pgRouting (for shortest path).
The road network I have does not include links for ferries, which is the preferable route in some cases, and is the shortest route predicted by Google/Bing etc. Our state (WA) has 10 ferries, so I can maybe manually connect the points to create links in QGIS and then potentially re-create the topology.
Is this a suggested approach for solving this problem? Is there any other easier or robust way to achieve this?
Another specific complication for my case is that the vehicles are electric, so the ferry section of the route, should not add distance to electric vehicle miles traveled. So, ideally the ferry section should have zero distance for all associated calculations. Is it possible to join the start and end points with a zero dist line? This would be analogous to shorting two points in an electric circuit.


texts:
I am looking to crop an ECW image based on latitude and longitude.
I have the following code to do that but it doesn't work.

<gdalwarp -t_srs EPSG:27700 -te 51.50386 51.50391 -0.4333941 -0.4333865 image.ecw imagecrop.tif>

I get the following error.

  Warning 1: -ts values have minx >= maxx. This will result in a horizontally flip
  ped image.
  Warning 1: -ts values have miny >= maxy. This will result in a vertically flippe
  d image.
  Creating output file that is -1297P x -1297L.
  ERROR 1: Attempt to create -1297x-1297 dataset is illegal,sizes must be larger t
  han zero.

How can I fix this?


texts: Spatial query in ArcGIS Javascript Api 3.1

texts: I have scenario in which I have to select features drawing a polygon and and get the attribute values in to a tables or lists.The features on the map consist of polygon, line,point at the same time so when he selects the features according to layer name I have to bind to separate tables.Is there any example tackling this scenario.


texts: Best method for expanding Multi* geometries to their constituent objects

texts:
I have a SpatiaLite table with a geometry column containing a mix of multi and single geometries, e.g. (Multipolygon, Polygon) or (MultiLineString, LineString). The application that will be reading this data unfortunately does not handle Multi* geometries well so I'll need to expand them using ST_Dump or something similar.
I assume it's better to move the resultant rows to a new table, correct me if that's not the optimal way. I want to make sure the new table has a unique primary key for each row (expanded and non-expanded geoms). There are other columns in the table besides the geometry column so I would like those to be duplicated for each Multi geometry that gets expanded. e.g.
Original table
+-----------+---------------+-------+--------------------------------+
| object_id |     city      | state |              geom              |
+-----------+---------------+-------+--------------------------------+
|         1 | Dallas        | TX    | MultiPolygon(((<coords>)))     |
|         2 | Los Angeles   | CA    | Polygon((0 0,0 1,1 1,1 0,0 0)) |
|         3 | New York City | NY    | Polygon((3 0,5 1,4 2,1 0,2 1)) |
+-----------+---------------+-------+--------------------------------+

Expanded table (Lets say the MultiPolygon is made up of 4 polygons)
+-----------+---------------+-------+--------------------------------+
| object_id |     city      | state |              geom              |
+-----------+---------------+-------+--------------------------------+
|         1 | Dallas        | TX    | Polygon((8 4,1 1,1 1,3 0,2 0)) |
|         2 | Dallas        | TX    | Polygon((1 6,1 1,1 1,1 0,1 0)) |
|         3 | Dallas        | TX    | Polygon((1 1,1 1,2 1,1 0,1 0)) |
|         4 | Dallas        | TX    | Polygon((1 3,1 7,1 1,1 9,1 0)) |
|         5 | Los Angeles   | CA    | Polygon((0 0,0 1,1 1,1 0,0 0)) |
|         6 | New York City | NY    | Polygon((3 0,5 1,4 2,1 0,2 1)) |
+-----------+---------------+-------+--------------------------------+

Could this be done and can I generate my own primary key? In this case the PK would be object_id.
I have a few methods in mind:

Use GeoPandas to manipulate the columns and add the new rows.

Extract the inner geoms (Shapely multipolygon.geoms)
Append the extracted geoms as new rows
Drop the current object_id column
Add a new object_id column
Write the GeoDataFrame back to SpatiaLite

Just use SQL(via GeoAlchemy) to perform the above steps from (1). I would imagine row_num() would be used to generate the new primary key column. Would this be too complex/not possible?

A concern with option (1) is the amount of memory that might be used if the table is massive since GeoPandas loads everything into memory.
Which of these options might be better, or is there another way of accomplishing this?


texts: How to immediately refresh ol.source.TileWMS

texts: How can I immediately update or refresh or redraw the ol.source.TileWMS ?


texts:
I'm trying to write a python script to load the Environment Agency (EA) Flood Zones WMS layers into QGIS (2.14). I've provided a link to the EA's data catalogue at the bottom of this question. The layers I'm hoping to use from the catalogue are "Flood Map for Planning (Rivers and Sea) - Flood Zone 2" and "Flood Map for Planning (Rivers and Sea) - Flood Zone 3".

I have used the code below, copied from a similar question (link at the bottom), but rlayer.isValid() is always returned as False no matter how I edit the url. I can take the url in the example code below and manually add the WMS layers with no issues.

Am I missing something in my url? Or is the issue somewhere else?

Extract from my code:

urlWithParams = 'url=http://www.geostore.com/OGC/OGCInterface?SESSIONID=-2049535737&INTERFACE=ENVIRONMENT&version=1.3.0&service=WMS&request=GetLegendGraphic&sld_version=1.1.0&layer=eainspire2011-wms-nat_floodzone2_inspire&format=image/png&STYLE=default'
rlayer = QgsRasterLayer(urlWithParams, 'EA Flood Zone 2', 'wms')
rlayer.isValid()
QgsMapLayerRegistry.instance().addMapLayer(rlayer)

EA Spatial Data Catalogue

Similar question - "WMS layer in a QGIS 2.4 stand alone Python script"

For anyone wondering the correct working links are below:

EA Flood Zone for Planning - Zone 2:
http://environment.data.gov.uk/ds/wms?INTERFACE%3DENVIRONMENT--86ec354f-d465-11e4-b09e-f0def148f590

EA Flood Zone for Planning - Zone 3:
http://environment.data.gov.uk/ds/wms?INTERFACE%3DENVIRONMENT--87446770-d465-11e4-b97a-f0def148f590