texts: New attribute creation using ogr and regex

texts:
I've got a folder full of dxf files with elevation values included as part of a long text attribute. I've written a regex expression to extract the elevation value, and would like to add it as a new attribute to an output Mapinfo file.

I'd prefer to complete this using command line and ogr if possible, as it's already doing most of the work to convert the files. Python alternatively. QGIS as a last resort.

Could this be completed using a SQL query in the ogr2ogr command?


texts: Converting Shapefile to DGN dos not display all geometry in MicroStation?

texts:
I am trying to convert a shapefile of FEMA floodplains from shapefile format to DGN format for use in MicroStation.
I'm using the "Export to CAD" tool in ArcGIS Desktop 10.4

The conversion appears to go well, and the resultant DGN draws just fine in ArcGIS. However, when loaded into MicroStation, much of the geometry does not appear. Its like in MicroStation the geometry does not exist.

It appears that any geometry with an attribute greater than 0 in the "Shape_Are" attribute field of the original shapefile gets ignored.

Can anyone tell my why these polygons are not being brought into MicroStation? Since the DGN draws correctly in ArcGIS, I assume that the geometry is written into the DGN but is not being recognized by MicroStation.
I am using:

MicroStation Select Series 2
ArcGIS Desktop Basic 10.4


texts: Batch Processing single band raster to create multi spectral images using arcpy.CompositeBands_management?

texts:
I am kind of new to python but I am trying to batch process single band raster images to create  multispectral imageries using arcpy Composite Bands. I have 24 single band rasters (Bands 1 - 3 ) all saved in a single folder which would create 8 composite images after its done. Fortunately the bands are in sequential order for each set, so I can input every three images without having to order them. So far I have:

arcpy.env.workspace = "C:\\Users\\Vinny\\Desktop\\Sentinel2_Scripts\\Unzipped\\Test_Out"

   direc = r"C:\\Users\\Vinny\\Desktop\\Sentinel2_Scripts\\Unzipped\\Test_Out"
out_raster = "C:\\Users\\Vinny\\Desktop\\Sentinel2_Scripts\\Unzipped\\Test_Out\\Composite_Bands"

i = 0
rasters = arcpy.ListRasters("*", "JP2")
for raster in rasters:
    i+=1
    arcpy.CompositeBands_management(raster, out_raster + "\\" + out)

Can you offer any suggestions or ideas?


texts:
I would like to do some word clouding as labels in QGIS from a series of tables in a database. Examples of word clouds can be found here: Word Cloud - Google Search

I have two tables in a database: locations and codes. The locations are cities that have geometry. The codes are the words I would like to use as labels. I have a count with the codes, which I hope to use as text font size.

Does anyone know if there is a way to label in QGIS from a related table. I have drawn a schema and example of what I would like to achieve.

Any thoughts are welcomed.


texts:
I posted a similar question as an answer on Stack Overflow, and it was recommended that I ask here instead.

I am developing an open-source library for converting between coordinate systems, and I've run into a wall.

I want to be able to perform a transformation on an MGRS point.  I want to be as accurate as possible, so I'd rather not convert to lat/long first, but that's the only way that I've been able to think about doing it.

Here are the circumstances:

I have a starting MGRS position
I see something in the distance and measure it's distance
I want to know the MGRS coordinate of said point

Is there any way of detecting whether I've crossed a zone boundary without converting to lat/long first?

It seems to me that MGRS zones (or UTM zones for that matter) are based on latitude/longitude.


texts: Serverless geospatial data type conversion (similar to ogr2ogr)

texts:
Is it possible to run ogr2ogr or Python/Fiona in a serverless environment where each feature gets processed in parallel and where you don't have to specify the properties of a virtual machine upfront?

I'm looking for a tool that can convert a variety of geospatial data types to something cloud GIS solutions such as AWS Athena, Google Biqguery and EarthEngine (in the future, current upload types are very limited) can work with.

An example I can think of is leveraging the power of Apache Beam but I have not seen any I/O connector capable of reading geospatial data types and create parallel tasks.


texts: Viewshed Analysis Memory Error in QGIS 2.18.15

texts:
I have been trying to run a Viewshed Analysis on a point layer containing approximately 174 points on a DEM compiled from the SRTM 30m Elevation data available at  http://dwtkns.com/srtm30m/

After I compiled a merged GeoTIFF from the required tiles I ran the Viewshed Analysis plugin on QGIS but after numerous attempts and changes in the input parameters I keep getting a MemoryError.

2018-01-02T12:20:46 1   Traceback (most recent call last):
          File "C:/Users/Frans_2/.qgis2/python/plugins\ViewshedAnalysis\viewshedanalysis.py", line 197, in run
            Algo)
          File "C:/Users/Frans_2/.qgis2/python/plugins\ViewshedAnalysis\doViewshed.py", line 409, in Viewshed
            mx_vis = numpy.zeros((full_window_size, full_window_size))
        MemoryError

Does anybody here have any idea why this is happening. I have updated the plugin to the latest available version (0.5.4) already. It seems other users have had similar errors but there is no clear answer to be found anywhere on how it was actually solved.

Project CRS is set to: WGS 84 EPSG:4326


texts:
Fairly similar to
Inside buffer with attributes in ArcGIS?
but in QGIS.

It would be convenient if the buffer layer inherited the attributes of the parent. Also I am doing inside buffers (-ve distance and dissolve).

If there is no such function, please let me know a workflow for joining the data.


texts: Deploying App to Portal for ArcGIS (ArcGIS Enterprise) from ArcGIS Experience Builder Developer Edition

texts:
I am new to ArcGIS Experience Builder Developer Edition but managed to create an experience with some theme customizations.
The next step would be now to deploy it to my Portal to make it accessible to my users.
I am trying to follow these instructions here: https://developers.arcgis.com/experience-builder/guide/deployment-topics/
However, it seems I am too much a beginner to get it work.
Is there any more detailed explanation available?
I am already lost on step one:

Unzip the downloaded experience and copy the folder to your web
server.

Where exactly should I copy the folder to?
I assume the web server is the server where Portal for ArcGIS installed but into which folder?
Also, are there more specific information on how to add and register the app?
The links provided refer to ArcGIS Online.


texts: Trouble Installing GriB extension to Geoserver 2.7.1.1

texts:
I am trying to add GriB reading capability to my Geoserver distribution (version 2.7.1.1, build date Jun 26, 2015). I am trying to use the plugin archive found at http://ares.boundlessgeo.com/geoserver/2.7.x/community-latest/

However, once I include all of the JAR files into .../geoserver/WEB-INF/lib/ and reboot the Geoserver instance, the Web Admin page returns a 503 error.

I have also tried copying all JARs save jettison-1.0.1 (to avoid overwriting the one in the main Geoserver distribution).

Would anyone have any insight why the GriB extension will not install properly, or which version of the extension I should be trying to install for GS 2.7.1.1?


texts:
I have a BMP file where I know the coords of the corners, and hence have a generated a World File according to ESRI's standards. It has the the right filename, but when I add the image in ArcMap, it complains there's no spatial information.

Things I've tried:

Turned on the "Use world file to define the coordinates of a raster" option
Different extensions for the world file - BMPW, BMW, BPW
Trying the help... the directions it gives are for a link file, which has to be imported manually for each raster (I'll have hundreds to these rasters)

So the question is: how do I get ArcGIS to use the world file?


texts: Select points with lowest elevations and same x,y coordinates

texts:
The problem is as following,

I have a point layer with x,y,z coordinates. Some of the points have same x,y but different z. I need to separate points without duplicates (x,y) and in case there are duplicates select only the ones with the lowest z.
How to separate them and have a layer with only 1 point per particular x,y?

I have ArcView 10.0 license.


texts: What causes ModelBuilder error occurred in script on this page MdDlgContent.htm?

texts:
I've been creating a model in Arc 10.0 for a few days now without problems. Then I added the Calculate Field Tool. I selected Make Variable --> From Parameter --> Expression. When I check off "Model Parameter" on the expression variable once it's been added to the model, save and close the model, then try to open it in Toolbox, I receive two popup errors, both saying "An error has occurred in the script on this page."

The first error says Error: Expected ')' in Line 79 Character 70 (Code 0):

The second says The value of the property 'IntForm' is null or undefined, not a Function object in Line 465 Character 1 (Code 0):

Both messages are referencing the following URL: file:///C:/Users/MYUSERNAME/AppData/Roaming/ESRI/Desktop10.0/ArcToolbox/Dlg/MdDlgContent.htm.

I located this article in ESRI support and followed the instructions, to no avail. (It must not be the exact same issue, as my message box doesn't reference Internet Explorer in the title). I also tried downloading SP5.

I opened the MdDlgContent.htm to take a look but that wasn't much help either because, despite some familiarity with Python and JavaScript, I am hardly a programmer.

Any help troubleshooting or resolving this issue would be greatly appreciated!!

I posted this same question in the ESRI Forums as well:
http://forums.arcgis.com/threads/79448-Model-Builder-An-error-has-occurred-in-the-script-on-this-page.-MdDlgContent.htm


texts: Problem processing Slovenian LiDAR DTM with GDAL (and QGIS): "too many stepY values" error

texts:
If I want to apply GDAL (version 2.2.3) on the LiDAR derived DTM (ASCII) from Slovenia (Link to Slovenian DTM below), I run into this error:

gdalinfo input.asc

ERROR 1: Ungridded dataset: At line 214, too many stepY values
gdalinfo failed - unable to open 'input.asc'.

The files can neither be loaded into QGIS (3.4.12).

I read here, that there is a XY swap in the files and sort can be used to get rid of the XY swap:

sort -k2 -n -k1 -t ';' input.asc -o output.xyz

However, after applying this to the .asc, I run into the following error:

gdalinfo output.xyz

ERROR 1: Ungridded dataset: At line 1002, too many stepY values
gdalinfo failed - unable to open 'output.xyz'.

What can I do to apply GDAL to the data?

The LiDAR can be accessed for free here: ARSO. Here is a direct download link to a sample ASCII file (D48GK projection): GK1_444_106.asc


texts: Issues while plotting the results of pgRouting shortest path

texts:
I am trying to plot the shortest path returned by pgRouting. However, I am having some issues. Sometimes it comes out a mess.

I have just joined the coordinates of the edges returned by the shortestpath query. However as mentioned in this link http://underdark.wordpress.com/2011/02/07/a-beginners-guide-to-pgrouting/ it says that it uses the necessary geometries. The path is plotted fine by QGIS. It is mentioned in this line

Therefore, we’ll join the results of the previous query with the network table containing the necessary geometries. Suggestions?

SELECT x1,y1,x2,y2
FROM shortest_path('
   SELECT id,
          source::integer,
          target::integer,
          cost::double precision
   FROM at_2po_4pgr',
   $start_node, $end_node, false, false)
INNER JOIN at_2po_4pgr
ON edge_id = at_2po_4pgr.id.

I used osm2po to generate the at_2po_4pgr table which contains the edges, their start and end vertices and other information.

I send the bunch of returned x1,y1,x2,y2 and then join these coordinates. Sometimes the path looks fine but sometimes it's a mess. The path looks fine when I use QGIS. But they have different implementation I guess.

Actually my question is, is it enough to send the start and end points of the edges returned by the query and then join these start and end points of each edge by straight line? Or I am supposed to use the geometry given by multilinestring which I have less idea about

I am trying to actually figure out what the "geom_way" column exactly represents and how to use that in my case. I created a network topology and persisted it to postgresql database using osm2po tool for a certain region. Now I am using pgrouting on it to find the shortest route. I am not sure how this column is used by qgis to display the route. It's just a bunch of number, how can I use it in my case?

I am actually using gpx format file to transfer the shortest path. How can I exploit this geometry and send that information in gpx so that it is displayed nicely.


texts: Changing legend item properties using ArcPy

texts:
I'm developing a Python script to created several map document (.mxd) files from .shp file and .xlsx file automatically. After setting symbology of layer and changing labels, I want to hide heading of legend by working on "show heading" option on "Legend Item Properties", but unfortunately I cannot found this property. Here is part of my code:

def CreateMXD(srcMXDDir, srcMXDName, srcSHPDir, srcLyrDir, srcStyleName, srcXlsxWB):
    srcLayer = arcpy.mapping.Layer(srcLyrDir)
    mxd = arcpy.mapping.MapDocument(srcMXDDir + "//" + srcMXDName + '.mxd')
    mxdDF0 = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    # mxdDF1 = arcpy.mapping.ListDataFrames(mxd2, "Layers")[1]
    mxdLayer = arcpy.mapping.Layer(srcSHPDir)
    arcpy.mapping.AddLayer(mxdDF0, mxdLayer, "TOP")
    addLayer = arcpy.mapping.ListLayers(mxd, "", mxdDF0)[0]

    fields = arcpy.ListFields(srcSHPDir)
    for field in fields:
        fieldType = field.type
        fieldName = field.name
        if fieldType in ["Double", "Integer", "Long"] and fieldName.startswith('edu'):
            arcpy.mapping.UpdateLayer(mxdDF0, addLayer, srcLayer, True)
            # print(addLayer.symbologyType)
            if addLayer.symbologyType == 'GRADUATED_COLORS':
                addLayer.symbology.valueField = fieldName
                labels = addLayer.symbology.classBreakLabels
                try:
                    addLayer.symbology.classBreakLabels = ["{} - {}".format(*r)
                                                                         for r in [
                                                                             [round(float(f), 2) for f in lab.split(" - ")]
                                                                             for lab in labels
                                                                         ]
                                                                         ]
                except:
                    print('Error in Symbology | %s' % fieldName)

This link has described this functionality in ArcMap but I want to do this change using ArcPy scripting. Below you can see the screen of my MXD file and mentioned option in ArcMap:


texts:
I was trying to measure a distance between a bunch of points and a polyline in QGIS and some people suggest me to use NNJoin Plugin. It worked smoothly but I crosschecked the results with  some manual measurement and in average the results are about 20 km off. Is there another way to do this with a more precise results?

edit: I try to measure the distance between volcanoes to its related subduction zone. I use WGS 84 Geographic coordinate system in my layers. I didnt do much with the NNJoin because there are only inputs for the two vector layer and for the output layer. The measurement that I used as reference is the measurement with measurement line tool.

actually I dont know for sure the unit of the NNJoin measurement results but I just asume its in hundreds kilos..


texts:
I use Oracle MapViewer to generate my mapping data and currently I'm using an Oracle database as my spatial back-end.  I would like to move my data over the PostGIS but I can't find any references if Oracle MapViewer actually supports PostGIS.

So my question is?  Is there out of the box support for PostGIS?  Would I have to write my own custom spatial data provider?


texts: How to solve "Cannot create topology" in ArcCatalog?

texts:
I have a FileGeodatabase in ArcGIS 10.0 with a feature dataset and five feature classes (3 polygon shapes, 1 line shape, 1 point shape).

I'm trying to create a topology to verify some vectors. Rightclicking Feature Dataset > New... > Topology.

I'm getting the following error message in ArcCatalog:

  ArcCatalog cannot create topology. The selected feature dataset does not contain any feature classes which can participate in a topology.

Browsing the internet for the error message, I found this tech article by ESRI. It tells the following:

  There are several reasons for a feature class to be excluded from participating in a topology:

  The feature class already participates in a topology or geometric network.
  The feature class is an annotation or dimension feature class.
  The feature class is registered as versioned.

I'm not quite sure what that means. There is no other topology or geometric network as far as I can see. What is an annotation or dimension feature class? The classes are all X;Y-coordinates only (2D). What is a versioned or registerd feature class?

I've tried checking all classes properties for something that sounds like that.

This answer to a similar question suggests to unregister the feature classes. But how?


texts: Is it possible to overlay ungeoreferenced raster images in the QGIS Desktop Browser?

texts:
I work for a company that produces Mining Investigation reports for property in the South West of England and I have spent the last 12 months digitising abandoned mine plans, geological survey, old ordnance survey maps etc onto a custom built GIS in QGIS. However, many of archive and analogue maps are extremely difficult to fit using a georeferencing tool, be it ARC or QGIS.

Therefore, it would be a lifesaver and an invaluable tool to be able to import and overlay ungeoreferenced raster images into the QGIS Desktop and be able to manipulate them over layers already in the browser, almost like a lightbox.


texts: QGIS print composer WGS'84 grid, distance between gridlines

texts:
I am trying to make an nautical chart in the QGIS print composer.

I have trouble getting the distance between my gridlines right. I need them to be a distance of whole minutes between. Now for some reason they are 1'2 apart. (I need them to be 003 46,00' - 003 47,00' - 003 48,00')

I tried with different CRS's and changing the interval in the grid settings, I can't get them right. What am I doing wrong?

Settings:


texts: Specifying Python code to be used on user form for QGIS?

texts:
Previously, on QGIS 2.18, I followed this well known tutorial to create a user form with some basic Python driven buttons.
https://nathanw.net/2011/09/05/qgis-tips-custom-feature-forms-with-python-logic/

The QGIS 3 layer properties interface no longer has the setting (at least not in the same place, so I can't find it) which tells QGIS to run a Python script on opening the form. This is discussed in the tutorial in the section headed "almost done".

How is this achieved?


texts: Cursor Update Python Script for Multiplying Several Fields against a Ratio (Real) in ArcGis

texts:
Disclaimer: Python Neophyte, Uber-reliant on GUI, working to wean...

I have a Python script for replacing multiple Field values with another table value that was graciously provided by the inimitable DKlassen and seen here:

import arcpy
tablename = "c:\temp\*NAMEOFYOURTABLE*"
fieldvalues = ['fieldname1','fieldname2','fieldname3'......'fieldnameX']

cur = arcpy.UpdateCursor(tablename)
for row in cur:
    sourcevalue = row.getValue('*name of field where source value coming from*')
    for f in fieldvalues:
        row.setValue(f, sourcevalue)
    cur.updateRow(row)

Question: After running it it seems to replace all the fieldvalues with f and defaults integers to 0 or 1 (since sourcevalues is a real number between 0-1)--how do I include an operator to multiple for f in fieldvalues by a ratio, or real/float number? Obviously the syntax of f*sourcevalue does not work.


texts:
I am trying to add a pie chart which summarizes data within each polygon based on two attributes fields ("P_Jet_2007" and "R_Jet_2007"). However, my script (see below) only produces an empty pie chart for each polgyon.

Please can anyone advise on what I am missing from my script?

uri = 'C:\PyQGis\P_Wood_Updated.shp'
lyr = iface.addVectorLayer(uri, '', 'ogr')

categories = [u'P_Jet_2007', u'R_Jet_2007']

colors = ['#3727fa', '#6810ff']
qcolors = []
for c in colors:
    qcolors.append(QColor(c))

canvas = iface.mapCanvas()

diagram = QgsPieDiagram()

ds = QgsDiagramSettings()

ds.font = QFont ("Helvetica", 12)
ds.transparency = 0
ds.categoryColors = qcolors
ds.categoryAttributes = categories
ds.size = QSizeF (100.0, 100.0)
ds.sizeType = 0
ds.labelPlacementMethod = 1
ds.scaleByårea = True
ds.minimumSize = 0
ds.BackgroundColor = QColor(0, 0, 0, 0)
ds.PenColor = QColor(100, 100, 100, 0)
ds.penWidth = 0

dr = QgsLinearlyInterpolatedDiagramRenderer()

dr.setLowerValue(0.0)
dr.setLowerSize(QSizeF(0.0, 0.0))
dr.setUpperValue(2000000)
dr.setUpperSize(QSizeF(1000, 1000))
dr.setClassificationField("Tot_Jet") # Attributed used to control size of pie chart

dr.setDiagram(diagram)

lyr.setDiagramRenderer(dr)

dls = QgsDiagramLayerSettings()
dls.dist = 0
dls.priority = 0
dls.xPosColumn = -1
dls.yPosColumn = -1
dls.placement = 0
lyr.setDiagramLayerSettings(dls)

label_settings = QgsPalLayerSettings()
lyr.setLabeling(QgsVectorLayerSimpleLabeling(label_settings))
lyr.setLabelsEnabled(True)
lyr.triggerRepaint

QgsProject.instance().addMapLayer(lyr)


texts: Can't have cubic resampleing with base map in QGIS 2.18.3

texts:
I want to display a raster with cubic resampling with base map, but it seems impossible in QGIS for some reason.
I first change (in the style tab of the raster) the resampling to cubic and I get the proper result:

Then, I add base map, for example google hybrid, and the raster changes to nearest neighbor automatically, even though in the style tab it is still on cubic. The image looks like that:

I've tried different base maps and I get the same results.. it changes back to nearest neighbor and I can't change it to cubic. Also, even if I remove the base map, I can't change it to cubic, I need to restart QGIS for that.

Is there a way around it so I can have both cubic and base map?


texts: QGIS 3.18. Export to GeoPDF - not showing label mask on contours?

texts: when exporting to GeoPDF, the label mask on contours does not work. This is very important to me. Is there a way to implement this?


texts: Importing points within a series of non-geolocated pdf maps and Excel data

texts:
I have a series of maps that have been drawn on by hand (points and areas) that have been given to me as PDFs, along with a series of interview PDFs explaining what those points and areas are (there are about 15 different categories). My job is to input all this into a map. I have created the overall map and was using the Freehand Georeferencer plugin to import the PDF maps, then resize and align them. Then I was placing the points by creating a new vector layer, one for each category, opening the attribute table and inputting the information from the interview pdfs. I would really like to be able to put this up on the web as well, or at least retain the capacity to do so.
I have three problems:

I have upgraded to 3.18 (Windows) from 3.12 and it appears that many of the points are no longer in their original locations. (For example, some points are in the middle of lakes, and I know they don't go there.)
I would like to simply input all the information into an Excel file, import that into QGIS 3.18 and then place the points by using the freehand georeferencer plugin.
The third problem is that it appears that the freehand georeferencer plugin now expects all pdfs to be georeferenced.

Can anyone give me any kind of guidance here? I know there are many potential complications, the least of which is my massive ignorance in Python. I'm in the first few lessons from "Learn Python 3 the Hard Way" and I can see that this will be a steep learning curve.
To make things more fun, this data is confidential so I have to be very general about any kind of specifics.


texts:
My task is to store the values in the database(postgreSQL) when the user inputs the lat-lon values in the form. and after submitting the form, the map should be displayed with markers at corresponding lat-lon values stored in DB. I have achieved this task using PHP,ajax & storing values in a text file. Now I have to do the same task using Geoserver+PostGis for which I am totally new. I have just installed Geoserver,postGIS,PostgreSQL in my system. I want your help to how to carry this task step by step. How to store these values in DB, how to retrieve them using geoserver.

Thanks in advance


texts: Knowing/converting AREA, PERIMETER units in QGIS?

texts:
I use Convexhull tool in QGIS and there are result value :

AREA = 0.0017494390 PERIMETER = 0.1725495172

I thought the result was strange. So, I took a rough estimate of the points through the services provided by other websites and picked out the results.

Other site result :
Area = 17220976.88 m^2 PERIMETER = 17030.657 m

Area value is a difference of 10^10
Perimeter value is a difference of 10^5..

I use "EPSG:4166 - Korean 1995" CRS and I also change project properties for units of distance, area. but, there are not any change.

So, I want to know what convex hull units are in QGIS and how to set other setting for ideal value.


texts:
As per - https://root676.github.io/installation.html - I am trying to install QNEAT3 via the OSGeo4W shell using python-qgis -m pip install matplotlib, however I am getting the following error message:

'python-qgis' is not recognized as an internal or external command, operable program or batch file.

I have had look and can't seem to find any info on this error, does anyone know why I am getting this or how I can fix it? I am running QGIS 3.4.1


texts: How to make my KML/KMZ accessible to the public

texts:
I created a basic web dashboard with MVC 4 that has a Google Map open on it. I have a KMZ file I made with ArcGIS desktop of a road in a Google Drive folder that i made publicly available so that i could take the URL of my file and implement that in MVC to overlay the KMZ file onto the Google Map, However when i run the process the dashboard opens up the Google Map on the correct spot however there is no shape file. The images attached show the code in MVC 4 that generates the Web dashboard, and the KMZ file sitting in the Google drive Folder. I am wondering if maybe there is not a problem with my script or if Google drive is the problem.

This is the URL to the KML: https://drive.google.com/file/d/0B4_TJdssPL1RRC13MTlzbjVPWHM/edit?usp=sharing

I am using MVC 4 and Google Javascript API v3


texts:
I've been struggling a lot lately to produce a map in R with the ggplot2 package. So far this is the map I managed to achieve (where each little plygon stands for a city) :

The data used for this map is here, the code used for this map is:

#require
library(ggplot2)
library(rgdal)
library(plyr)
library(dplyr)

#open df
data = read.table("datafr_comm.txt", header=T, sep="\t", quote="", dec=".", colClasses=c(rep("factor",)))
data$score <- as.numeric(as.character(data$score))

#open shapefile
mapa <- readOGR(dsn="France",layer="codes_postaux_region")
mapa <- spTransform(mapa, CRS("+proj=longlat +datum=WGS84"))
#prepare for merging dataframe/shapefile by ID
mapa@data$id <- rownames(mapa@data)
mapa@data   <- join(mapa@data, data, by="ID")
mapa.df     <- fortify(mapa)
mapa.df     <- join(mapa.df,mapa@data, by="id")
#merge dataframe/shapefile
plotDatafr <- left_join(mapa.df, data)

#plot map
ggplot() +
geom_polygon(data = plotDatafr, aes(x = long, y = lat, group = group), fill = "white", colour = NA) +
geom_polygon(data = plotDatafr, aes(x=long, y = lat, group = group, fill=item), colour = NA) +
scale_fill_manual(values = c("#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"), name = "") +
coord_map()

Now, I can't figure out a way to fill automatically the white zones on the map by calculating an average of the surrounding points. I tried many things including the stat_density_2d or stat_bin2d functions, but no way to achieve what I really want.


texts: WFS in OpenLayers4 returns data but doesn't display any

texts:
I can't get this WFS to work.
I tried to stick to the OL-WFS examples - no luck.

let vectorLayer = new ol.source.Vector({
  format: new ol.format.GML(),
  url: function(extent) {
    return 'https://geodienste.hamburg.de/HH_WFS_Statistik_Stadtteile_Wahlergebnisse' +
      '?version=1.1.0&request=GetFeature&typename=Statistik_Stadtteile_Wahlergebnisse:Buergerschaftswahl_15.02.2015_-_Wahlbeteiligung_in_Prozent;
  },
  strategy: ol.loadingstrategy.bbox
});

this.featureLayer = vectorLayer;

return new ol.layer.Vector({
  source: vectorLayer,
  style: new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: 'rgba(0, 0, 255, 1.0)',
      width: 2
    })
  })
});

This query returns 1.5mb of data - but nothing gets displayed.
I also tried the same with this format:

new ol.format.WFS({gmlFormat: ol.format.GML2()})

As well as new ol.format.GML 2 & 3 - no luck.
And I tried it with this line of code at the end of the url:

'&bbox=' + extent.join(',') + ',EPSG:6.9:25832'

If I add this (as seems to be recommended in the OL examples) the 1.5mb of data is not returned. It should be a very simple thing to include this data in an OL map!


texts: Many-Many Relationship in GDAL Pixel Functions

texts:
I am using pixel functions (in Python) in VRT derived data sets. As far as I know it is only possible to use One-One or Many-One relationships within a pixel function. In other words, the in_ar includes one raster or multiple, whereas out_ar is flat, i.e. it should be a 2D array.

Is it somehow possible to save multiple output bands or is there a workaround?

As an example: I want to sort values of the 3D input array along the third dimension and then store each layer/level to a new output bands.


texts:
There are two rasters A and B. A is smaller but within B and since A has more up-to-date pixel values, I was able to merge A and B with "Mosaic raster layers" tool in QGIS which replaces B pixel values with A. However there is a problem, that A does have no value pixels (which are dummy high numbers e.g. -3.40282e+038), when the merge happens it overwrites the values of B with high noData values which I don't want it to happen as B has values for those missing points within A. Is there any way it could be done in QGIS? It should be possible with gdal_calc.py, but I don't know how to set up the rules.

I have tried the following so far, but i hope there is a one-step way of doing it:

Clipped larger B to the extent of A (otherwise gdal_calc.py doesn't work for different dimensions)
gdal_calc.py -A A.tif --NoDataValue=-3.40282e+038 -B B.tif --outfile=result.tif --calc="where(A==-3.40282e+038,B,A)"

After doing these steps i was planning to do mosaic of result.tif and B.tif, however result.tif is not correct, so I must be doing the where rule in gdal_calc.py incorrectly.

Rather than doing a lot of steps, there might be one tool which specifies "do not merge NoData of A" etc, but at the moment the merge product is coming back with -3.40282e+038, whereas original A raster didn't show up this value in symbology.


texts: Using Create Fishnet and Project in ModelBuilder?

texts:
When I try to run the model below, the Create Fishnet tool creates a huge grid that totally ignores the extent input (BLOCK_Project1). It's nowhere near the extent input either.

If I click on the Project tool and run it independently, then do the same with the Create Fishnet tool once the Project tool is shaded, it works and gives me this.

I just can't seem to run it all at once and get the results I'm looking for. Whenever I do, I get this for the fishnet.

Any idea what's causing this issue?


texts: Earth Engine how to avoid timeout for buffer calculation?

texts:
I have a function that takes a feature and creates a buffer around it.

createBuffer = function(feature){

  var AOI = feature.geometry().buffer(3000);

  return ee.Feature(AOI).copyProperties(feature);
}

When I run this function on a collection that contains polygons that have  coastlines, the fractal nature of the coastline makes the buffering take a while, the result is delivered to another function to do imageCollection calculations based on the new buffered feature.

In some cases (though none that I can replicate here) it leads to a The operation timed out.

What can be done to avoid these timeouts?


texts:
I am trying to update attributes in a table using Calculate Field in Model Builder. I have used the following bit of code but its not updating the field.

    expression:
    Change(!Zone!)

    def Change(Zone):
      if (Zone == "NULL"):
      return "Non UK"


texts: What type is this PostGIS data and how can I get lat, long from it?

texts:
I have an R dataframe, with location stored in a column. Location values look like this:
0101000020E6100000797A00E6B5593E40D5CE30B5A5394940
I was told that it is PostGIS data.
The goal is to visualize all those locations in R. My plan is to use ggmap, and for that I need to have lat, long values.
I have spun up a PSQL DB with PostGIS extension to try and poke at this data, but I have no clue what type it is, and treating it as varchar seems like a dead end. I have looked into PostGIS documentation on data types, but haven't been able to find anything that resembles my data. My initial guess was that my data is of point type, but it doesn't seem to be correct, since points are represented as pairs of numbers (for 2d), not as 51-chars-long strings.
Now I am at a loss. What type is this data and how can I get lat, long from it?


texts: Using TIFF file in Intersect tool of ArcMap

texts:
I am attempting to do some geoprocessing with a CSV file and a TIFF file (trying to use the Intersect tool), but am experiencing some trouble with it. I can't seem to be able to input the TIFF file as a layer for the tool.
Do I have to convert the file into another format?
I have never had any formal education in ArcGIS and am teaching myself.


texts: Unique and Range value symbology in Arc Desktop 10.1

texts:
I am facing a issue in which i have to give symbology depending on unique value of one attribute  and range value of another attribute.
That is symbology on gender and age.
like

symbol1 - Male(gender) and 21-30(age)
symbol2 - Male(gender) and 30-50(age)
symbol3 - Male(gender) and 50 -100(age)
symbol4 - Female(gender) and 21-30(age)
symbol5 - Female(gender) and 30-50(age)

Thank you in advance.


texts: QGIS identify tool not working on selected layer

texts:
In QGIS For some reason, the identify tool doesn't work on the layer I have selected in the TOC. Instead, it just selects whichever the largest overlapping layer is. For example, I want to identify a point that's within a larger polygon, it will always select the polygon even if the point layer is highlighted. So in order to use identify to check the point, I have to deselect every layer that's overlapping the point layer. It's a bit annoying since I work with several land classification layers that overlap one another.
This only started a few months ago, any way to.


texts: What software for direct editing of a DBF file?

texts:
I am currently processing GPS log data from my holiday and need to add information to the points. I remember I used OpenOffice and later LibreOffice to directly edit the files but currently I am really upset because I always get a tiny DBF + huge DBT file instead of only one DBF. QGIS is then not able to load the shapefile correctly... :-(

Tried various software but no success. Even old OpenOffice 2.4 which I currently downloaded and ran using Wine does the same problem. Gnumeric is not able to save the table, Sdbf probably does not support UTF-8...

My reason for manual attribute editing is that I cannot comfortably edit them in QGIS - it was much easier for example to mark beginning and ending point and than copy the atributes in spreadsheet...

Have no idea why it is not possible yet. Tried win machine with newest LibreOffice and unfortunately the same problem.

any ideas?

thanks


texts:
I have two columns one that has the BeginLigeSpanVersion and one that has the end LifeSpanVersion, both kind of noted into Date Time format. I said kind of because usually in the unfiltered unaltered state it looks like this in string format:
{ "DisplayValue": "2020-02-05T16:37:15.4534428+01:00", "Value": "2020-02-05T16:37:15.4534428+01:00" }
With the magic of this command it is converted to date time :
CASE WHEN  "beginLifespanVersion" Is NOT NULL
THEN
to_datetime(string_to_array("beginLifespanVersion"  ,'"')[7])
ELSE
NULL
END

Now it took me a while to realise that the difference between the two dates is two months...kinda. I have posted a picture how it suppose to look like below. How do I do how do I add two months to this?


texts: I have created an ArcGIS enterprise geodatabase. There are i5, i6,... i30 tables in the sde schema of my geodatabase. What do these tables do?


texts:
I see from the ExportWebMap specification that rotation can be set in the JSON object sent to the sever.

Does anyone know if this can be accessed via the ESRI JavaScript API (v3.9)?

I do not need to display the map with a rotation, just honor it in a printout

Thanks,
Joe


texts: Does ModelBuilder need to use variable to define projection and then re-project?

texts:
I am trying to build a Model in ModelBuilder to create two different process raster datasets:

One set is to define a raster dataset to set a  WGS 84 projection; and then
the second set is to take the define projection raster datasets to reprojection to NAD UTM 83 to a new folder.

Do I need to create a variable to copy from a define projection raster set and then
Reprojection raster dataset to NAD UTM 83?

What is the correct variable tool should I use to perform to connect them?

BTW, I have to do this with a lot of raster datasets and do the same thing.


texts: QGIS 3.16 GeoTIFF Not Producing World File

texts:
I am using the function to export map to Image, I have the "Append georeference information (embedded or via world file)" ticked but I do not get a world file and when I drag the TIFF in to a GIS, it comes in at 0,0,0. i.e. there does not appear to be any embedded georeference.
How can I ensure that both the TIFF and the world file are created?


texts:
I have installed QGIS 2.0.1, great thing, but the OpenLayers plugin doesn't work correctly.  When I try to load two different layers I get the answer "could not draw Google satellite because:..."

Also there are problems with loading one map. The correct display is only possible in the OpenLayers projection EPSG:3852, when I try to change the project projection the display gets wrong. This is depending on the zoom factor of the display. At the same time the OpenLayers collapses.

I set the project properties on "on the fly"...


texts:
I have been trying to export a number of scanned paper maps onto my database, and have been just adding them as raster layers as part of a larger group layer. I georeferenced them, and none of them overlap (they are smaller towns around the county). When I create a mosaic dataset and add the rasters to it, the file ends up being around 40 GB. This is odd, because I also used the same process with about 500 maps of the same file size (which slightly overlapped and filled an entire area) and the final mosaic dataset only took up about 5 GB.

Is there a way to make the file size smaller or store the rasters more efficiently?


texts: Is there an se_toolkit version for ArcSDE 10?

texts: I have the one for 9.3 downloaded--is there a version for ArcSDE 10?  So far my 9.3 version of se_toolkit is working against 10 geodatabases, but I've only tried a couple of the tools(sdeping and sdestats).


texts: How to rename feature classes based on field values within the attribute table?

texts:
My goal is to iterate through several feature classes and rename each feature class based on a field value within the feature class's attribute table.  These feature classes are a single polygon feature, thus they only have a single row in the attribute table.  I was thinking I could use the Rename GP tool to do this so I've got that part figured out.  Next I can use the SearchCursor to find the field value that I want to use to rename the feature class with.

So I'm struggling here in two places.  First, I'm unsure how to iterate through each of the feature classes to get their current names and then apply that to the first parameter of the arcpy.Rename_management tool.  I know I can get the name of the feature class using arcpy.Describe but I'm stuck after that.  Second, I'm struggling with how to assign each of the field values from the attribute table to a variable to be used for the second parameter of the Rename tool.  I can easily print the name using rows = arcpy.SearchCursor(fc, "", "", "name")  for row in rows: print row.name but I'm struggling with how to assign its value to a variable.

Let me know if I need to explain this any better.


texts:
I have to assign the feature type of a layer.
Does anyone know an function, or way to access this information.

I have noticed that the message of the function readOGR (package rgdal) contains an entry like Feature type: wkbMultiPolygon.


texts: Finding specific pixel values in ArcGIS using Python

texts:
I have a raster file containing some integer values.

I need to find all pixels whose values are not defined in list like this one:

l = [1,2,4,8,16,32,64,128].

Is there any easy way to solve that problem without converting raster to points and editing attribute table?


texts: How do i use gdal2tiles for creating tilecache?

texts:
I want to create tiled cache from raster file, so i downloaded gdal2tiles.py.
Now i don't get what (and how) i must do with that python script file.


texts:
I wanted to clarify with all of you whether I am setting up Cesium correctly in wampserver using Windows 7 64 bit.
The Reason on why I want use wampserver is because I want to access data from the sql database(localhost/phpmyadmin/) and php script.
I wanted to clarify whether I am doing correctly before I proceed to starting my project.

I am going to list down the steps below.
I am using Windows 7 64 bit

Download wmapserver and install on C
I download git and git clone git clone https://github.com/AnalyticalGraphicsInc/cesium.git to wamp64/www/
I type npm install in node.js command prompt
And then I type npm run release
I navigate to the Services utility and set wampapache64,wampmariadb64,wampmysqld64 to Automatic under StartUp Type
If I want to run the cesium file, I type the link http://localhost/cesium/Apps/HelloWorld.html in browser to start without typing node server.js

Am I doing the correct method to set up Cesium in wampserver? Also did I make any mistake? If there is another correct way to set up Cesium in wampserver?

My Requirement
I want the Cesium website to access the data from the database sql.
I want the Cesium website to access the php script.

Can you help me ?


texts: QGIS sum up polygon features inside grid

texts:
I have a grid with population points inside each of the individual cells. I've used the "join attributes by location" tool to join each population point to their respective cell. Now I'd like to sum up all my population points inside each grid cell. I've used the field calculator using the following expression:

sum("id")

"id" being the id of each cell in the grid. But now, I have a sum of population for each row.
I'd just like to have one line with the cell id and the total population inside.


texts: Least cost path analysis : a lot of sources and destinations

texts:
I'm working on a 15*10 km area, splatted in a 50*50 m cell grid. I want to compute the time needed to travel from a single starting point to each cell centroid, combining driving and walking time.

I placed a point every 30 m along my trail/roads. Using pgrouting,I already have the least time cost from my single startpoint to each of these points.

I computed the shortest path of each centroid to the nearest trail/road point. Assuming the same point-centroid walking speed (4 km/h), I get a time for each. I had first used the "Walking time" extension on those segments, but some times weren't realistic (canyons).

But the time costs limits are raw on the raster: the nearest point on the road/trail isn't necessarily the fastest.

Beside that, this approach doesn't take into account the slope, I tried r.walk and r.drain, but I can't figure out how to use it with all this : there are over 4 000 roads/trails nodes and 27 000 centroids !

How can I get the least-cost path from any of the centroid to the fastest trail/road point ?

I'm using Qgis and Grass.


texts: Using two Search and Update cursor pairs on the same feature class

texts:
I've been working on this code for about 3 days and all of my time has been spent trying to understand why I can only execute the deleteRow() method on the first pair of cursors.  The idea is that a cursor will iterate through the rows of a .dbf table and remove any rows whose attributes values in a dictionary created by the SearchCursor.  Here is the applicable code:

with arcpy.da.SearchCursor(ParcelJoin, ["Name"]) as cursor:
    lookupTbl = {}
    for row in cursor:
        lookupTbl.update({row[0]:row[0]})

with arcpy.da.UpdateCursor(Foreclosure, ["Parcel"]) as cursor:
    for row in cursor:
        if row[0] in lookupTbl:
            cursor.deleteRow()
        else:
            continue
    print('Joined features have been removed!')

with arcpy.da.SearchCursor(AddressJoin, ["FullAddr"]) as cursor:
    for row in cursor:
        lookupTbl.update({row[0]:row[0]})

with arcpy.da.UpdateCursor(Foreclosure, ["Address"]) as cursor:
    for row in cursor:
        if row[0] in lookupTbl:
            cursor.deleteRow()
        else:
            continue
    print('Joined features have been removed!')

Foreclosure is a .dbf of property sales, ParcelJoin is a join between a parcel layer and the Foreclosure table based on ParcelID and AddressJoin is the same feature class joined on address.

The first iteration of the cursor pairing executes successfully and removes the relevant rows.  However, the second pairing iterates properly (I've checked by creating a counter to see how many times the code iterates) but does not execute the deleteRow().  I've looked all over the web but can't find any examples of code that use two cursor pairs in this way.

Finally figured it out, there were trailing whitespaces in the Foreclosure table that needed to be stripped. Code with list comprehension works beautifully!


texts: How do I reverse geocode Census Block group from 15 character FIPS Code to LATLONG?

texts:
I know that the block group is a polygon but the Census Gazetteer files provide a lat/long for census tracts here - http://www.census.gov/geo/maps-data/data/gazetteer2010.html

Where can I get this for the census block group level ?

The FCC provides an API that converts LATLONG to CENSUS BLOCK GROUP but Id like the reverse [https://www.fcc.gov/developers/census-block-conversions-api]- maybe Im being vacuous ?


texts: Why ArcGIS Runtime SDK for .NET Toolkit is made as open source?

texts:
Is there any specific reason for making only ArcGIS Runtime SDK for .NET Toolkit an open source project where as the actual ArcGIS Runtime SDK for .NET is proprietary?

Though this question sounds weird, I just want to understand whether to deliver more useful features based on the community feedback?


texts: Reading attribute popup information from cluster presented points in the map generated via qgis2web?

texts:
I have a problem with a web map generated via qgis2web plugin using clusters. I have a batch point data and many of these are on the same location. Because of that I have used cluster for representing the locations. The problem is that when I click on the cluster and got the tree with points at the end and after that click on specific point on of these points doesent showing the popup window (if I have 4 points, one of these is always bad, the popup windows does not appear).  Marked point on the screenshot is that point with popup error. I have tried to change cluster circle radius in css folder but I have not found the code part for it.

How do I solve this problem?


texts: Why do TIF loaded into PostGIS lose their colormap?

texts:
I have successfully loaded a TIF image (using raster2pgsql) into my PostGIS 2.0.1 database. I can load the raster image into QGIS and it displays in the correct space. However, it has lost all colour and appears as a greyscale image. (Layer Properties > Style > Single band gray, colormap = grayscale)

If I load the same image directly from file into QGIS it displays colours correctly.(Layer Properties > Style > Single band gray, colormap = colormap)

If I convert the TIF to JPG and load the JPG to PostGIS and then display in QGIS the colours render correctly. (Layer Properties > Style > Three band color)

I can save the colourmap from an image and apply it to the greyscale images but the colours aren't quite right.

What do I need to do to get the TIF images to retain their colourmap when I load to PostGIS?

Windows 7 32-bit
PostgreSQL 9.1.4
PostGIS 2.0.1
QGIS 1.8.0
raster2pgsql -s 27700 -I -C -M *.tif -F -t 200x200 gis_rasters.osvmd > vmd.sql

Thanks in advance

Ross


texts:
I have two street data sources, and as you can see, they are not perfectly aligning. However, I need data from both of them, so I wonder if there is a way to make them align perfectly with each other. Also,  I want to fix A due to its accuracy (assume one is called A and the other B), and make B follow the shape of it.

I am using ArcGIS Desktop 10.3.1 Advanced License.


texts: Getting "repeat values" individually to calculate them separately

texts:
I'm working on smart form of maintenance activities by XLSForm and trying to calculate the weightage of each question to evaluate the priorities of taking action based on the final score of the survey ether it is High or Low.

The engineer want to have weightages in each inspected assets like Transformers, RMU etc.

I designed the questions of these assets to be in the "repeat".

It might be two or more transformers inside the substations so the repeat will increase.

How can I get repeat values individually to calculate them separately?

1

]2


texts: Read negative elevation data from NASA png tiles

texts:
I query NASA JPL service using query string. Response is a unsigned 16 bit per pixel png image and the values are elevation values in meters. Same is described in the capabilty xml of the layer as can be seen below:

<Layer queryable="0">
      <Name>worldwind_dem</Name>
      <Title>SRTM derived global elevation, 3 arc-second</Title>
      <Abstract>
    A global elevation model, prepared from the 3 arc-second SRTM dataset by filling some of the problem areas. Prepared by the NASA Learning Technologies.
    The default style is scaled to 8 bit, non-linear.
    It is possible to request the elevation data in meters by the short_int tyle and requesting PNG format. The resulting PNG file will be a unsigned 16 bit per pixel image. The values are then the elevation in meters.
    Values are signed 16 bit integers, but PNG will present them as unsigned, leading to a few areas with very large values (65000+)
    For elevation values in feet, request PNG format with the style feet_short_int.
      </Abstract>

When I read the elevation using following code

var pixel = bmp.GetPixel(rowIndex,
                         colIndex);
var value1 = BitConverter.ToInt16(new byte[] { 0, pixel.B },
                                  0);

I read the positive elevation value correctly but unable to read negative value. The query string tile mentioned above covers all the case, it has data for +ve, zero and -ve values.


texts: Field calculator to calculate if field contains certain text

texts: I have a field (WS_Type in the screenshot) that I want to use to calculate values based on the text within the field. The number of occurrences of each WS_Type in the field is not important, but the presence of each different type in the field is. Basically, I need to create an if/then statement that says, if this field contains "NRW" AND "Diversion", then it will get a score of X. If this field contains "NRW" OR "Spring" it will get a score of Y, etc. I'm not sure how to set up this if/then statement.


texts: How to densify LineString vertices in Shapely/GeoPandas

texts:
I'm looking for an equivalent of densify tool in ArcGIS, but with GeoPandas or Shapely. I have a line that would like to be densified with a specific distance; how to add vertices of a line using Python?
Suppose I have a line that has 2 vertices (WKT)
LINESTRING (0 0, 9 0). how do I divide the line into an equal distance of 1  such that it returns the same line, but with instead of 2 vertices, it has 10 vertices? the result I expect:
LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0, 6 0, 7 0, 8 0, 9 0). I find this post which is similar to what I want, but it uses ArcGIS instead of GeoPandas and Shapely.


texts:
I cannot get to open my Python scripts saved in default directory: ....\QGIS\QGIS3\profiles\default\processing\scripts in QGIS 3.0 Processing Toolbox. There is no folder called Scripts like in 2.18.

Anybody knows how to open my Python script in Processing Toolbox of 3.0?


texts: Is there a way to group regions on a world map for common functionality?

texts:
I'm new to trying to use mapping, and not finding the kind of documentation I need.  I'm looking at OpenLayers, but open to other solutions.

I have a client that has sales reps with territories worldwide.  Those territories are not country exclusive, there may be more than one rep in a country with territories divided between states/provinces.  I need the map color-coded so the territories are highlighted by sales rep.  OpenLayers at least, looks like it might be the way to go because zooming in on a world map more sharply defines the regions within a country  (zooming in on the US, shows the individual states).

As it stands, I don't know if OpenLayers will support splitting a country on an un-zoomed view, or how to color code any of the territories (many of the example pages for OpenLayers don't seem to work).


texts: How can I use QGIS 3.8 to work with features from the new Geological Map of South America (2019)?

texts:
The Columbian Geological Service has recently uploaded a fantastic new geology map for all of South America:

There are multiple file formats available for download. You can work with the map in ArcGIS formats (mxd, gdb), print (pdf, tiff), and web viewing options. However, I can't get the files to load in QGIS.

QGIS should open the GDB file, but the format isn't recognized. I can load the raw shapefiles from that folder, but all the style and detailed color information is lost. So how can I get these map features into QGIS 3.8?


texts: PostGIS - Error when converting shapefile into dataset using openJUMP

texts: I am new to this field. I want to save my shapefile as a dataset in PostGIS using openJUMP, but I'm getting an error. Can anybody suggest a solution for this?


texts:
I'm using a .csv detailing types of renewable energy in the UK. I want to display them as different colours for the types of energy (e.g. wind, solar, gas etc.) but also a different size symbol based on the energy output.

All the information is in the .csv file and I am able to perform these two styles independently, but I can't work out how to do this together.  I'm new to GIS so just getting the hang of things!


texts:
I thought that I could open KML layer like a normal vector layer but it doesn't work (it says "invalid datasource"). How can I open a KML layer in QGIS? Or how can I open it anyway, I never worked with kml-data before...
I'm working with QGIS v. 2.0.1, the OS is Windows.


texts: QGIS: Dissolve after raster-to-polygon not working correctly, only dissolving a portion of the vector layer

texts:
I am working with QGIS and am trying to convert a raster layer into a vector layer.  I have converted it to a polygon from raster (image below), and then wish to dissolve this polygon into a single vector shape. This has worked for the majority of the shape, but at some longitude there is an error, and it only dissolves vertically but not horizontally (2nd image below).

Has anyone had this error before, or can anyone recommend a different approach to get a vector layer from a raster?

I have used the Dissolve function, and also tried the GDAL/Grass dissolve functions with no success.


texts: XMLHttpRequest cannot load error When we have GeoServer and WebServer on same domain

texts:
I have Geoserver and webserver on same domain (Windows XP with IIS). But when i tried to get WMS feature info, i got this error in chrome. If Geoserver and webserver on same orgin, no need of proxy i think. so please suggest me to overcome this issue. Thanks in advance

XMLHttpRequest cannot load http://localhost:8080/geoserver/cite/wms?LAYERS=cite%3Aplym_p1_examplec…e%2Fpng&INFO_FORMAT=application%2Fvnd.ogc.gml&SRS=EPSG%3A27700&X=596&Y=160. Origin http://localhost:3890 is not allowed by Access-Control-Allow-Origin.


texts: How to display join query results on a sublayer in CartoDB.js application?

texts:
I'm working with my first CartoDB.js application.

I created a visualization on the app using

cartodb.createVis('map', visjson_url, opts);

The app has two layers, base map and data layer.
the data layer has five sublayers. Each sublayer represents data in an independent table (i.e., sublayer1 represents table1, etc.).

I run queries using this statement:

data_layer.getSubLayer(sublayer_index).setSQL(query_str);

The problem is that I get an error message saying "Bad Request" if query_str involves join operations between some of the tables (e.g., spatial join between table1 and table3). If query_str is related to only one table (without joins), setSQL() above works fine.

How can I display join results using setSQL() that involves join operation between two or more sublayers?


texts: Is there a way to track edits which have been made to a feature class in ArcGIS Pro? I.e. a field being auto populated which would give information in regards to what kind of edit has been made, e.g. "[Feature] location has altered" or "[Feature] [Field name] has been altered"


texts: Finding the closest points to a line with only one point from each subset of points


texts: Using output of unsupervised cluster analysis to train supervised classification for vegetation mapping?

texts:
I'm using Google Earth Engine to classify mangrove forests from Landsat imagery to calculate extent. Let's say I perform an unsupervised K-means cluster analysis on a Landsat scene to identify general patterns and themes in the spectral bands.

After grouping the clusters together into several relevant classes (e.g. mangrove, other_vegetation, soil, urban, water), would it then be possible to use this classified image for training a supervised classification algorithm such as random forest?

The idea is to use a stratified random sampling technique (~5000 points) to assign the training points within the classes obtained through the cluster analysis.

Can anyone see a flaw in this method and would it increase the accuracy of the final classification?

TLDR - is it possible (or advisable) to use the output from an unsupervised classification technique to train a supervised classification algorithm?


texts: Bulk import geotagged raster images in attribute table

texts:
I have more more than 700 points in my shapefile with name of images in attribute table. I want to add a raster field to the shapefile so that i can add all images of those points in respective field.

Adding manually is tedious.

Can you suggest a better way for me to import those geotagged raster into attribute table?


texts: Extract data from Geopackage .gpkg using lat long bounds

texts:
I have a .gpkg.tar.gz data file whose size is 3.8 GiB. Uncompressed to .gpkg, its size is 8.5 GiB.

It is spanning a large geographical area. I want to extract data from this for a small area, given by lat-long bounds.

It is fine if there are features in/crossing this area but extending beyond it - I don't want to bother with chopping/clipping features. But I want to exclude features that are entirely outside of the given bounds.

With .pbf (Openstreetmap) data I have a way of doing this, example:

osmconvert india-latest.osm.pbf -b=78.55834,12.54397,80.3558,13.71654 --complete-ways -o=small.pbf

Is there a similar command / process we can do with .gpkg?

I understand it is originally SQLite; and I've read that one doesn't need to load the entire file into memory to work with that. So some specific queries - table names, column names might be there? This file is too big for me to load up in QGIS etc at the moment to be able to investigate its contents.

Ref: http://www.geopackage.org/index.html


texts: Making polygon shapefile from unsorted point shapefile?

texts:
I have a point Shapefile with unsorted feature ID's (FID). I want to convert the same point file into polygon Shapefile by using python gdal/ogr. I got the help from link1 and link2. I'm able to make polygon shapefile successfully but because of my unsorted point Shapefile the output of polygon shapefile is irregular (not as expected). Expected polygon shapefile is ring/polygon/box around points .

.
Is there any way to sort point shapes as right to left or left to right before making polygon Shapefile?


texts: Can GDAL be used to load raster data into Postgres/Postgis geodatabase?

texts:
I would like to know if GDAL can be used to load raster data to a postgres/postgis geodatabase

I had read that there was problems with the GDAL postgis raster driver?

I was wondering what particular GDAL commands would be used to load the raster data to the geodatabase

I am not sure if GDAL would allow me to load each raster image at a time or require a mosaic to be produced first.

Are there any parameters that can be used/set to improve loading speed?


texts: Collapsing layers in ArcMap Table of Contents using ArcObjects?

texts:
I'm stuck on how to collapse all layers in my ArcMap projects (9.3). I've tried a few bits of sample code, but they all seem to fail for some reason or other.

Some of the layers might be grouped, and some are imagery or WMS layers if that makes a difference.

Here is some code I have tinkered with, but does not work:

Dim pMxdoc As IMxDocument
Set pMxdoc = ThisDocument
Dim pFeatLayer As IFeatureLayer
Dim pEditLayers As IEditLayers
Set pEditLayers = pEditor
Dim pMap As IMap
Set pMap = pMxdoc.FocusMap
Dim pLayer As ILayer

'Collapse any expanded grouplayers
Dim pGroupLayer As IGroupLayer
For i = 0 To pMxdoc.FocusMap.LayerCount - 1
Set pLayer = pMxdoc.FocusMap.Layer(i)
'Check if the reached layer is a group
If TypeOf pLayer Is IGroupLayer Then
    Set pGroupLayer = pLayer
    if pGroupLayer.Expanded = true Then
        pGroupLayer.Expanded = false
    End if
End If
Next i
pMxdoc.UpdateContents
pMxdoc.ActiveView.Refresh


texts:
Firstly, apologies as I'm a newbie and have only been studying GIS for a few months.
I am digitising field boundaries and have a field within a field, and am wondering what the best way of representing this in a single layer is, so that the area of the internal field is subtracted from the area of the main field.  I have read an explanation about donut parcels but am not sure if this is the right solution for what I want to do? I am using ArcMap.  Thanks in advance for your suggestions.


texts:
Is there a possibility in QGIS to populate a ValueMap-like form widget with predefined data and additionally offer a possibility to fill in data (that might be missing) manually? Like an editable combobox.
EDIT:

Did you mean adding missing value while setting up your form widget or while editing attribute?

When I enter my data it would be nice to select predefined values, but if by chance a not predefined value appears in field-work it would be nice to record it as well.
A possible solution would be, when an entry like "other" is selected another textbox will be enabled to type the "not predefined value".
But the drawback of this solution would be, to have an extra field instead of using the same.
QGIS-Version: 3.16.5 LTR


texts: Finding centroid for bounding box (bbox) in R

texts:
I'm trying to find the center of a bounding box created using the function bbox, the main reason for that is that get_googlemap() requires a center for it to plot the map.
I'm trying to plot a set of CSV data with Lat and Long, with the google background. I'm already able to acquire the maps from google using the API's key.
I've tried to work with the average of my lat and long columns, but the results weren't so good.
The coordinates were acquired in datum WGS84.


texts:
In OpenLayers 4 changelogs, it explains:

  Unlike the switch from v2.x to v3.x, which marked a complete rewrite
  of the library with an entirely new API, major version increments now
  simply mean that users should pay attention to the 'Breaking changes'
  section of the upgrade notes.

But "semanting versioning" says:

  MAJOR version increases when you make incompatible API changes

So, I don't understand: if the OL API is not changed, why OL major version did it increase to 4 ?
Can someone explain better?


texts: Export multiple Images from collection to Google Drive not working with Google Earth Engine Python API 0.1.238

texts:
I am struggling with Google Earth Engine Python API version 0.1.238 to export images from the Landsat EVI collection (from 2002 to 2020, but i started with a small sample) to my Google Drive, which are essential to my research.
At first i wanted to add these images directly in QGis, but apparently the version of the Earth Engine plugin of QGis is outdated for such task (i.e. export + impossible to loop the Map.addLayer command).
So I finally built a python script prompting no errors, but my Google Drive folder is still empty. Here is what gives the task.status() command after the task:

{'state': 'READY', 'description': 'myExportImageTask', 'creation_timestamp_ms': 1602830108152, 'update_timestamp_ms': 1602830108152, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': '2QGWJCVBPUOOCTNFIIDYLKRQ', 'name': 'projects/earthengine-legacy/operations/2QGWJCVBPUOOCTNFIIDYLKRQ'}

How could i make this script work ?
import ee
ee.Initialize()

# USDA NAIP ImageCollection
collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_32DAY_EVI')

# create a roi
polys = ee.Geometry.Polygon(
        [[[37.5, 25.5],
          [40.5, 25.5],
          [40.5, 28.5],
          [37.5, 28.5]]])

# filter the ImageCollection using the roi
naip = collection.filterBounds(polys)

# filter dates
naip_dates = naip.filterDate('2002-04-01', '2002-08-07')
mosaic = naip_dates.mosaic()

# print out the number of images in the ImageCollection
count = naip_dates.size().getInfo()
print("Count: ", count)

# export the ImageCollection to Google Drive
img_lst = naip_dates.toList(100)

for i in range(0, count):
    image = ee.Image(img_lst.get(i))
    name = image.get('system:index').getInfo()
    task = ee.batch.Export.image.toDrive(image=image,
                                         region=polys,
                                         folder='32DAY_EVI',
                                         scale=30,
                                         crs='EPSG:3857')
    task.start()
    print(task.status())

I also tested such simple script to test my Google Drive and still nothing in my folder (not even the creation of the folder itself):
import ee
ee.Initialize()

landsat = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_123032_20140515').select(['B4', 'B3', 'B2'])

geometry = ee.Geometry.Rectangle([116.2621, 39.8412, 116.4849, 40.01236])

task_config = {
    'scale': 30,
    'region': geometry,
    'folder':'GEEtest'
}

task = ee.batch.Export.image.toDrive(landsat, 'TEST_todrive', **task_config)
task.start()

I am wondering if it is not an issue with my Credentials. I did the earthengine authenticate step without any issue, but is there something else to do to link GEE to GDrive ?


texts: How do I make a constant slope elevation model in arcMap?

texts:
I am currently working with ArcMap to develop an elevation model of a quarry showing how the terrain at the quarry site would look post-closure. One of the requirements is that the slopes placed along the benches should be 30 degrees exactly. I have an elevation model of the site and a shapefile showing all the quarry leases. One of the assumptions is that the boundary elevation shall remain the same in the future and all of the rock within the boundary would be extracted. The final base elevation is assumed to be 100m for the time-being; i.e. elevation of the pit surface post closure.

To start, I extracted the elevation points of the boundaries where the slopes shall be put in; using "Extract by Mask" and "Sample" to tabulate the values and then import them into a new template. I was wondering if it would be possible to create the points at the base of the slope using Excel; i.e. design some sort of an algorithm that will determine the XY coordinates of a point inside the boundary such that:

The distance between a point on the boundary and the corresponding point inside the boundary is (Boundary Elevation - 100)/tan 30deg.
The point inside the boundary is such that it's on a line perpendicular to the line where the boundary point is situated.

However I'm pretty sure this would take quite a lot of time and would not be the optimal way of going about the problem. As the final base elevation is subject to change, it would be better to have a solution that doesn't require much effort.

What would be the optimal way of creating the elevation model?

I have gone through some of the forums in ESRI as well as stackexchange but haven't been able to understand the methods suggested, particularly those involving python as I don't have any experience working with Python. So I was wondering if it will be possible to work on it without using python or if it does require the use of python, could anyone explain to me how the script can be written?

The method suggested for this question for instance: Creating constant slope DEM using ArcGIS Desktop?


texts: Shape Normalization Error while importing shapefiles in cartoDB

texts:
I had succesfully installed CartoDB in Ubuntu 14.04 LTS VM and could import CSV and GeoJSON without any error. When I try to load shapefiles, the following error pops out.

Error Code : 3009

Error : SHP Normalization Error

My environment for python is Version 2.7.6 and I have latest GDAL installed which is of version 2.1.3.

I have tried to export the uploaded GeoJSON file and then try to import in CartoDB to prevent any encoding/projection issue but still no success.


texts: Is there any documentation on how to change the user interface of Potree?

texts:
Has anyone found or created any information regarding the different user interfaces of Potree?

I have created LiDAR portals using the LAStools laspublish tool to create Potree sites, but cannot see how I could adapt the GUI.

Here are two examples:

http://potree.org/demo/potree_1.3/showcase/tern_auscover_chowilla.html
http://potree.org/showcase/lion_head_simone_garagnani.html


texts: hub_distance in python standalone script

texts:
I need to load geometry from csv and then use it to calculate hub distance. I write the following script, but it doesn't work. What am I doing wrong?

 sys.path.append('/usr/share/qgis/python/plugins/')
 from processing.core.Processing import Processing
 from processing.tools import *
 qgs = QgsApplication([], True)

 QgsApplication.setPrefixPath("/usr", True)

 Processing.initialize()
 QgsApplication.initQgis()
 layer = QgsVectorLayer('from_excel.csv', 'from_excel', "delimitedtext")

 layer.setCrs(QgsCoordinateReferenceSystem(32365, QgsCoordinateReferenceSystem.EpsgCrsId))

layer2 = QgsVectorLayer('/home/raven/work/UpWork/Inputs/Cenefas.shp', 'Cenefas', 'ogr')
layer2.setCrs(QgsCoordinateReferenceSystem(32365, QgsCoordinateReferenceSystem.EpsgCrsId))

general.runalg("qgis:distancetonearesthub", layer, layer2, 'test', 0, 0, 'HubOutput')

QgsApplication.exitQgis()

Then I've bunch or errors:

Warning: Not all input layers use the same CRS.
Fields of width 255 truncated to 254 'NoneType' object has no attribute 'boundingBox'
Segmentation fault (core dumped)


texts: How to Center a Split Colour Ramp on Zero in Arcgis?

texts:
I have a number of different rasters with different ranges. I want to display them with defined intervals of 10. Values range from -200 to 400 varying for each raster. I want to symbolize them with the same colour ramp. Colour ramp below.

The issue is that the transition of the split colours does not always fall on zero. Is there any way to make the transition fall on zero without having to do it manually every time?


texts: Creating loop for Lines data (spatial data)?


texts: Exporting an image in Google Earth Engine image alignment

texts:
When I export an image in Google Earth Engine the resulting raster does not match the original. When I import the export back into Google Earth Engine the is an spatial misalignment between the original (orange in picture below) and the exported version (green in the picture below). The direction as size of the misalignment varies across the image. A solution appears to be to set the scale of the image below its native resolution, in this instance 1m as opposed to 30m; however this will hugely increase the size of the image. Am I missing something in the export call

Export.image.toDrive({
  image: landsat,
  description: 'imageExport',
  scale: 30,
  crs: 'EPSG:3857',
  maxPixels: 800000000000
});


texts: Clip only returning top polygon in ArcPy?

texts:
I'm running the command below on a data set with many stacked polygons.  This is a county standard so I can't change the way they manage their data.  Each type of polygon can be uniquely identified with one of the fields (TaxLot, Block, ROW, etc.)

arcpy.Clip_analysis('parcels', 'clip_sect', parcel_clip3_shp)

It works, but it only returns the top most polygon. However, running it as a geoprocess within ArcMap returns all the polygons.  I built the same geoprocessing tool in ModelBuilder and it returned all polygons.  I exported the code to *.py and it only returned the top polygon.

I was thinking I could create feature layers for each type with a selection query and then run a function to clip/add but that seems a bit brute force and not very clean.

Any thoughts?

I am using ArcGIS 10.2.1, Python 2.7.3 and Idle 2.7.3.


texts: GRASS is crashing every time I am running this pyqgis script. Any suggestions?

texts:
I have the following script that it keeps on crashing. How can I find a solution to this issue?

import pandas as pd
import processing
from PyQt5.QtCore import *
from PyQt5.QtGui import *

resol=0.0007

uri = 'file:///'+FolderPath.csv?delimiter=,&crs=epsg:4326&xField=longitude&yField=latitude' csv_file = QgsVectorLayer(uri, FolderPath, 'delimitedtext') QgsProject.instance().addMapLayer(csv_file)

parameters_Red = {
    'input' : csv_file,
    'npoints' : 1000,
    'power' : 2,
    'column' : 'Value',
    'rcolumn' : '128:000:000',
    'GRASS_REGION_PARAMETER' : 'Xcoord1, Xcoord2, Ycoord1,Ycoord2 [EPSG:4326]',
    'GRASS_REGION_CELLSIZE_PARAMETER' : resol,
    'GRASS_RASTER_FORMAT_OPT' : '',
    'GRASS_RASTER_FORMAT_META' : '',
    'GRASS_SNAP_TOLERANCE_PARAMETER' : 1,
    'GRASS_MIN_AREA_PARAMETER' : 0.1,
    'error':'TEMPORARY_OUTPUT',
    'output' : 'FolderPath/GRASS_OUTPUT1.tif'}

processing.runAndLoadResults('grass7:v.surf.idw', parameters_Red)

`


texts: Google satellite xyz tile layer in QGIS has lower resolution than the basemap in Google Earth?

texts:
As far as I can tell, the Google satellite imagery loaded as an xyz tile in QGIS is the same image that is used in Google Earth, but at a lower resolution. The lower image below is from QGIS using the xyz tile while the top image is from Google Earth. Both are screenshots.

I used "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}" for the xyz tile and I am using QGIS 3.8, but I also have this issue when using older versions. Is there a way to improve the xyz tile resolution?

Edit 1: the CRS for the QGIS image is WGS 84/UTM zone 17N (EPSG,32617)
Edit 2: they are just simple windows snips, but do not look any different when exported using 600 DPI


texts: Identifying and deleting overlapping or near lines with similar running?

texts:
I have shapefiles with lines that shows surface water runoff. I need to clear and simplify these shapefiles by condition that if some lines are near by other and has similar running, this lines has to be deleted.

I'm using Python 2.7, ArcMap 10.2 and Windows 7.

I have tried two solutions like:

Solution 1:
Iterate through lines using ESRI arcpy module, make buffer with 10 meters distance on each line separately and then use Select By Location with Completely Contains condition (without the original line from which it was buffer made) and delete this lines. Then iterate next line. This works fine, but Select By Location was TOO slow, and because I have 244 shapefiles that needs to be cleaned and together contains 48 milions of lines, so it approx. 200 000 line per shapefile and I need script to do it quicky.

Solution 2:
Like in previous solution, I use iteration through lines and separately use buffer, but then I want to avoid to use Select By Location, so I use Clip the lines with buffer and then compare the lengths of unclipped lines and clipped lines. If the lengths was the same - so it means that lines completely belong into buffer - I sum the IDs of this lines to list and delete from original shapefile (no from the clipped, it was only temp). But in this solution, there was Clip function so damn slow, so it is not usable for me.

Now I want to leave using arcpy module, because it not suitable for using in such a large data. I want to use GDAL/OGR to reconsturct a procedure of one of previous solution, but I'm really loosing in OGR algorithms. All that I founf out till now is load shapefile and iterate through lines, but I can't figureed out how to make buffer, identify nearby lines with similar running and delete them.

Pictures showing iteration on one line:

And pictures showing final result:


texts: Render QgsRasterLayer without adding to legend

texts:
I am trying to render a number of related QgsRasterLayer objects that are related to eachother (they are slippy-type tiles). Ideally I would like one item in the legend that can have its position changed (i.e. move above or below other legend items) and add custom actions/properties dialog to, but have the actual rendering of the tiles be done by QgsRasterLayer.

My first approach was to attempt to subclass QgsPluginLayer and override the draw() method to call QgsRasterLayer.draw() for each of my loaded rasters. This does nothing, even if layers are added to the QgsMapLayerRegistry.

def draw(self, renderContext):
    for rasterlayer in self.loadedlayers.values():
       rasterlayer.draw(renderContext)

My second approach was to create a group in the legend that is somehow hidden. I can add the layer to the layer tree using

mygroup = QgsProject.instance().layerTreeRoot().addGroup("mygroup")
for rasterlayer in self.loadedlayers.values():
    mygroup.addMapLayer(rasterlayer)

but it doesn't seem possible to hide this group from the user. From Showing only some layers in QGIS legend? it doesn't seem like there is an alternative. I feel as though there is an easy solution to this that I'm not seeing.

EDIT: I tried creating a QgsMapLayerRenderer object that renders the layer to an image using a QgsMapRenderer, and while this occasionally works, it often produces output that appears to have elements of the current screen in the output image (maybe having something to do with using a QgsMapRenderer() while one is already invoked on the current thread?). In my QgsPluginLayer I have:

def createMapRenderer(self, context):
    return MultiRasterRenderer(self, context)

And my QgsMapLayerRenderer is this (includes some debugging output):

class MultiRasterRenderer(QgsMapLayerRenderer):

    def __init__(self, layer, context):
        QgsMapLayerRenderer.__init__(self, layer.id())
        self.layer = layer
        self.context = context

    def render(self):
        fout = open("/Users/dewey/Desktop/outputdump.txt", "w")
        rendererContext = self.context
        if len(self.layer.loadedlayers) > 0:
            img = QImage(QSize(self.context.painter().device().width(),
                                    self.context.painter().device().height()),
                              QImage.Format_ARGB32_Premultiplied)
            painter = QPainter()
            painter.begin(img)
            painter.setRenderHint(QPainter.Antialiasing)

            render = QgsMapRenderer()
            render.setLayerSet([layer.id() for layer in self.layer.loadedlayers.values()])
            render.setProjectionsEnabled(True)
            render.setDestinationCrs(rendererContext.coordinateTransform().destCRS())

            fout.write("Context extent: " + rendererContext.extent().toString() + "\n")

            if render.setExtent(rendererContext.extent()):
                render.setOutputSize(img.size(), img.logicalDpiX())
                render.render(painter)
                painter.end()

                img.save("/Users/dewey/Desktop/output.jpg")

                rendererContext.painter().drawImage(0, 0, img)

                fout.write("Render extent: " + render.extent().toString() + "\n")
            else:
                fout.write("render.setExtent() returned false\n")
        else:
            fout.write("No tiles to load, skipping rendering")

        fout.close()
        return True


texts: Why does Germany look equally as large as Afghanistan on Google Maps?

texts:
I snipped Germany and Afghanistan from the same Google Maps view:

They look like they are about the same size. Here is the full screenshot:

The problem: Afghanistan has an area of  while Germany has an area of  (according to Wikipedia).

Now, I know that Germany appears larger because higher latitudes are stretched more when drawn on 2D. But I calculated the difference between the perimeter on latitude 34° (Afghanistan) and latitude 50° (Germany).

Given that the Earth's radius at the equator is 1, the relationship between latitude and the radius on that latitude should follow the circle function

 gives us y = 0.926
 gives us y = 0.831

So the radius is 0.926 in Afghanistan and 0.831 in Germany.

This is only a difference of 11% caused by the higher latitude. According to Wikipedia, however, Afghanistan is 83% larger than Germany.

So where does the rest of Germany's scaling on Google Maps come from?

Please excuse my incorrect naming of geographic terms, I have absolutely no experience in geography.


texts:
In the PostGIS documentation it says that there are two steps to creating a spatial table with SQL:

  Create a normal non-spatial table.
  Add a spatial column to the table using the OpenGIS "AddGeometryColumn" function.

If I followed the examples, I would create a table called terrain_points like this:

CREATE TABLE terrain_points (
  ogc_fid serial NOT NULL,
  elevation double precision,
);

SELECT AddGeometryColumn('terrain_points', 'wkb_geometry', 3725, 'POINT', 3 );

Alternatively, if I look at existing tables in pgAdmin III, it seems like I could create the same table like this:

CREATE TABLE terrain_points
(
  ogc_fid serial NOT NULL,
  wkb_geometry geometry,
  elevation double precision,
  CONSTRAINT terrain_points_pk PRIMARY KEY (ogc_fid),
  CONSTRAINT enforce_dims_wkb_geometry CHECK (st_ndims(wkb_geometry) = 3),
  CONSTRAINT enforce_geotype_wkb_geometry CHECK (geometrytype(wkb_geometry) = 'POINT'::text OR wkb_geometry IS NULL),
  CONSTRAINT enforce_srid_wkb_geometry CHECK (st_srid(wkb_geometry) = 3725)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE terrain_points OWNER TO postgres;

-- Index: terrain_points_geom_idx

-- DROP INDEX terrain_points_geom_idx;

CREATE INDEX terrain_points_geom_idx
  ON terrain_points
  USING gist
  (wkb_geometry);

Do these two methods produce the same result? Is the version based on pgAdmin III simply more verbose, and doing things that AddGeometryColumn would do by default?


texts:
A simple to ask but difficult to answer question:

I want to create a control like the one in the Leaflet layers example:

http://leafletjs.com/examples/layers-control.html

Notice that the control is a small button, until pressed. When pressed, it expands, and shows all options. I can't find anywhere in the documentation on how this was implemented, or how I can copy this behaviour. Does anyone have any idea?


texts: NetCDF time aggregation using R stars package

texts:
I have a monthly temperature netCDF dataset with aprox 20 years and I'd like to calculate a monthly climatology (average of all Jan, Feb, Mar...).

I know I can do this with stackapply from package raster. But I'm trying to learn the package stars. I tried using the aggregate function but I'm not able to figure out how to set the by argument correctly to apply the mean function to each month.

This is what my dataset looks like:

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
 Tavg_monthly_UT_Brazil_v2_19800101_20151231.nc
 Min.   :18.11
 1st Qu.:25.00
 Median :26.36
 Mean   :25.85
 3rd Qu.:26.85
 Max.   :29.83
 NA's   :56236
dimension(s):
     from  to offset delta  refsys point                    values
x       1 168 -74.25  0.25      NA    NA                      NULL [x]
y       1 162   6.25 -0.25      NA    NA                      NULL [y]
time    1 432     NA    NA POSIXct    NA 1980-01-16,...,2015-12-16

Is aggregate the correct function to use? I also tried group_by but that does not work for stars objects.


texts: How do you set up XY Event Layer Parameters on ArcGIS Script Tool? (asking for user input)

texts:
I've made the XY Event Layer work through python script, but can't get it to work correctly as a python script tool on arc, asking for user input. Does anyone have any suggestions?

This is the script I am using on arc:

# Import arcpy module
import arcpy, sys, traceback
from arcpy import env

# Script arguments
Input_Table = arcpy.GetParameterAsText(0)
feat_output = arcpy.GetParameterAsText(1)

# Local variables:
XY_Event_Layer = Input_Table
Longitude_and_Latitude_Points = XY_Event_Layer
sr = arcpy.SpatialReference(104145)
feat_output = "long_lat.shp"

# Process: Make XY Event Layer
arcpy.MakeXYEventLayer_management(Input_Table, "longitude", "latitude", XY_Event_Layer, sr)

# Process: Copy Features
arcpy.CopyFeatures_management(XY_Event_Layer, Longitude_and_Latitude_Points, "", "0", "0", "0")

This is my parameter set up of the script tool properties:

input Table   -  Table View,
x coor     -     Field,
y coor     -     Field,
output feat.  -  Feature Layer


texts:
I converted a csv file to projected SpatialPointsDataFrame in R, input it to bkde2D, converted the result to raster, saved the raster to a TIFF, and input that raster into the GME Isopleth tool.

Here is my R code:

    data = read.csv("_PolarBearData.csv", na.strings = "NA", header = T)
    oldproj = "+proj=longlat +ellps=WGS84 +datum=WGS84"
    myproj = "+proj=stere +lat_0=90 +lon_0=-55 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"
    coords = cbind(data$longitude, data$latitude)
    spdata = SpatialPointsDataFrame(coords, data, proj4string=CRS(oldproj))
    # ensure it's projected
    spdata = spTransform(spdata, CRS(myproj))

    # KDE parameters
    xy = cbind(coordinates(spdata)[,1], coordinates(spdata[,2])
    bandwidth = 34390
    grid_x = 130
    grid_y = 288

    # Run KDE and save to raster TIFF
    datakernel = bkde2D(xy, bandwidth, gridsize=c(grid_x, grid_y))
    ras = raster(list(x=datakernel$x1, y=datakernel$x2, z=datakernel$fhat))
    writeraster(ras, paste("mydata", '.tif', sep=''), "GTiff", overwrite=T)

Next I ran this GME code to produce 95% volume contours:

isopleth(in="I:\Data\mydata.tif",
         out="I:\Data\Contours.gdb!mylines",
         quantiles=c(.95),
         poly="I:\Data\Contours.gdb!mypoly");

Here is an image of my output: Bad GME Isopleth

Black-to-white raster (background) represents a  KDE generated in R with "bkde2D" from the KernSmooth package. The yellow points were the input to the KDE. The red line shows the 95% isopleth generated in GME. GME is stretching the contours south when the input is an untampered TIFF from R, but otherwise, if the raster is from Arc for instance, GME produces correct contours. Why does this happen?

UPDATE:
I should also mention that I added the raster from R directly into ArcMap, and it had an unknown CRS even though I specified it in R (at least, for the input). I tried using Define Projection to fix this, but GME still produces the same output. The raster itself has coordinates that make sense (it matches the point clusters in the image), but GME is not producing good isopleths from it.


texts:
I'm a python newbie.

I have a Geotiff file. I need to find coordinates in the Geotiff to find the position of a thermometer within the geotiff (the thermometer has the coordinates -115.47083484975207, 35.75837483423 e.g.).

I have these geotiff information:

PROJCS["NAD83 / UTM zone 12N",
GEOGCS["NAD83",DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4269"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-111],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","26912"]]

How do I transform the geotiff, so that I get coordinates in lat and lon? I have already tried this

But I don't quite know if this leads me to a solution. I just want to know how to find to position of the thermometer in the geotiff.


texts: Using the ArcGIS Runtime SDK for C# is there a method to Hit Test for a polygon on the graphics layer?

texts:
I am very new to using ArcGIS runtime SDK for C#.
I have loaded a polygon from a shapefile.
I wish to query if a lat/long point is inside or outside the polygon?

How can I do this?

Thanks in advance!


texts: Has anyone integrated Nearmap imagery into a CMV application?  Have you utilized a Tiled Map Service (TMS), WMS or something else?  Have you run into any problems with a need for credentials especially if incorporated into a public application?


texts: How to set access to GWC urls in GeoServer?

texts:
I published a service in GeoServer successfully and now I want to set user access to my demo url (http://myip:8080/geoserver/gwc/demo).
Everybody can see my demo url without any login and grant and this is my issue. How can I limit access to the url with authentication?


texts: QGIS: Intersect two polygons, one line layer as Output

texts: I have two polygons that intersect, now I want to calculate the length of the outline of the new polygon, but only the outline it shares with one of the input polygons. Is there a way to enter two input polygon layers and receive one, specifically requested output line layer?


texts:
I have a bit of trouble visualizing my data.
I have a dataset with 60 minutes of seismic signals and i determined the signal source areas. My goal now is to find out how the results change using different window sizes.

Example: The dataset is divided into 5 minute windows and the signal source areas are calculated for each window (i.e. each window is treated as its own dataset). I then overlap the results of each window to see which signal areas are detected using this window size. The localization count shows how often the areas overlap.

I did this for every window size and the end result looks as follows:

Now i want to merge these subplots into one plot, showing which signal source areas appear in which window size. So basically it should show the areas that the window sizes have in common and which are different. And you should be able to recognize which window size(s) each source area belongs to.
I don't know how to show this in a clear, well-arranged way.

Maybe i'm using a wrong approach, i feel like this could be handled very well with some spatial technique. I did some research but i couldn't find anything that could be applicable.

The localization counts can be omitted, it's just about the area shapes now.


texts: Fill/interpolate gaps in raster area WITHOUT extending the outer borders

texts:
I have a raster file that is related to a river bed in a curvy river section. Inside this raster file are gaps with no data (e.g. small islands, measurement gaps, etc...).
I want to fill those gaps inside the raster. However, I couldn' figure a way or find a function that fills nodata gaps without extending the outer border (see image). I tried interpolating the dataset but this returns interpolated values extending the original border.
Is there a function that I oversaw and anyone here experienced a similar issue with raster data to be filled/interpolated?
Thanks in advance


texts: Generating MapService tiles cache from machine other than server with ArcGIS Server?

texts:
I have a MapService published on ArcGIS Server v10.3. I would like to generate the tiles cache for this service. For an unknown reason ArcGIS Server CPUs are getting overloaded to 99%-100% after a few day of processing.

When that happen, all the other services published to ArcGIS Server are not responding very well until I kill the tool that generate the tiles. Then, I can restart the tool and continue to generate the tiles for another few days before the CPUs get overloaded again.

I would like to know if there is a way to generate the tiles for a MapService from another machine, like a virtual machine with an ArcGIS Desktop license and a writing access to the server drive disk.

The generate time is not really important. I would like to prevent the tiles generation process to have a negative impact on all the services.


texts: Is the erratum about polygon point order in GeoJSON correct?

texts:
The GeoJSON specification (RFC 7946) states that for polygons:

  A linear ring MUST follow the right-hand rule with respect to the area it bounds, i.e., exterior rings are counterclockwise, and holes are clockwise.

However, this erratum published in July 2017, states:

  It seems like the clockwise/counterclockwise descriptions are the opposite of the right-hand rule. Walking an exterior ring in a counterclockwise direction would have the exterior of the ring to the right of the observer.

This is a pretty significant issue, that does not appear to have been resolved (judging from the "reported" status of the erratum).

Is this correct? What exactly is the "right-hand rule" anyway?


texts:
I've got a list of some hundreds of latitude/longitude pairs close to each other (within a country). I want to create an image with a map, and circles overlaid around these points. Basically, what I need is:

an image of a map
a way of translating latitude/longitude pairs to pixels

Where can I obtain a fitting map (e.g. roadmap or satellite map) with the parameters needed to plot the points? I suppose there is some FOSS tool that does most of this already?


texts: I want to submit a paper for some research that I have done. What are the top conferences and journals in this field?


texts: PostGIS returning a field one level higher in the Json

texts:
I have the following SQL on a PostGIS Database ($properties contains a list of fields, including osm_id)

        SELECT row_to_json(fc)
            FROM ( SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
            FROM (SELECT 'Feature' As type,

            ST_AsGeoJSON(ST_Intersection(lg.way, bbox.geom))::json As geometry,
                    (
                       SELECT row_to_json(t)
                         FROM
                         (
                              SELECT $properties
                         ) t
                    )
             As properties
          FROM
              planet_osm_polygon As lg,
                (
                    SELECT ST_Transform ( ST_SetSRID(ST_MakeBox2d( ST_MakePoint( $topLon,$topLat) , ST_MakePoint( $botLon , $botLat) ), $srid) , 900913 ) as geom
                ) as bbox

                WHERE "building" is NOT NULL AND ST_Intersects(lg.way, bbox.geom)
                    ) As f
                ) fc;

that returns GeoJson in this format

"type": "FeatureCollection",
"features": [{
    "type": "Feature",
    "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
            [
            ],
     "properties": {
        "osm_id": 278325855,
        "height": "80.4"
    }
},

I need to return it in THIS format, note the extra id field at the properties and geometry level. The value to be returned is the same as the osm_id in the properties.

"type": "FeatureCollection",
"features": [{
    "type": "Feature",
    "id": 3698396,
    "properties": {
        "height": 54
    },
    "geometry": {
        "type": "Polygon",
        "coordinates":
            [
            ]
        ]
    }
},

How can I do this in a Query, or would it be more efficient to parse the Json and add it later?


texts:
I am attempting to convert a shapefile to WKT using ArcGIS 10.4. I have seen the question/answer posted here, Converting geometry to WKT using ArcPy?, which is helpful. However, when I follow the code and write it to the text file I create, the spatial reference is incorrect. The geometric shape is identical to the original shapefile, but the WKT maps to a different continent. I don't have experience writing to WKT, so I could be making an amateurish mistake that is causing my spatial reference to be off.

import arcpy

# Create new empty text file
outFile = open(r"my_path.txt", "w")

# My shapefile
fc = r"my_path.shp"

# Convert to WKT and write to file
cursor = arcpy.da.SearchCursor(fc, ["SHAPEFILE@WKT"])
for row in cursor:
    print(row[0])
    outFile.write((row[0]) + "\n")
del row
del cursor
outFile.close()

The projection that I want my WKT to be in is NAD_1983_UTM_Zone_13N.