texts:
I'm not sure where the right place to raise this is, so:

The PROJ strings for EPSG codes 4151 and 4283 are identical:

+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

which is a hassle because one is Switzerland and one is Australia. This causes an issue when writing geospatial data to file from R with sf, raster or sp and then trying to open it in other software like QGIS. 4151 is usually autodetected as the file CRS, presumably because the EPSG code shows up first in a search.

With PROJ 4.9.1 on Windows, I thought I'd found a way around this by including +init=EPSG:4283 in the string:

+init=EPSG:4283 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

However, with PROJ 5.2.0 on Ubuntu 18.04 LTS, this string is considered invalid and an error, 'no arguments in initialisation list' is returned. This is encountered both with sp::CRS() and sf::st_crs() and persists even when +init is used alone, e.g.

sf::st_crs('+init=EPSG:4283')

Validity appears to be case-sensitive - If I run

sf::st_crs('+init=epsg:4283')

no error is returned in Ubuntu. However, if I do that under Windows, I get inconsistent results:

> st_crs('+init=EPSG:4283')
Coordinate Reference System:
  EPSG: 4283
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
> st_crs('+init=epsg:4283')
Coordinate Reference System:
  EPSG: 4283
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
 > st_crs('+init=EPSG:4283 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
Coordinate Reference System:
  EPSG: 4283
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
> st_crs('+init=epsg:4283 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
Coordinate Reference System:
  No EPSG code
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

Note the missing EPSG code.

Is my best option as an R package developer to force the init tag to lower-case if a Linux-based OS is detected? I don't think I have an easy way to detect the PROJ version in use.


texts: Replacing NetCDF file under the hood in GeoServer

texts:
I'm receiving NetCDF files hourly, and would like to update them and then refresh them (as the time steps and values change).
Within the API, I'm not sure which endpoint would be useful after reading the documentation: https://docs.geoserver.org/latest/en/user/rest/index.html
Looking through the documentation I'm unsure of the proper endpoint to call.
My end goal is the replacement of the file on the filesystem followed by a call to refresh the layer with the existing parameters.


texts: Help using offline QGIS API documentation

texts:
I have downloaded the QGIS documentation for offline usage, am I correct in thinking that the qgis.qch file is the one I need to open to view the documentation, using Qt Assistant?

I have a Qt installer, where it lists a bunch of packages to install but Qt Assistant isn't listed.

(i did check if this has already been asked/answered - as far as I can tell it hasn't)


texts:
I used standart license initialization code in C# standalone program. It was working perfectly on 10.1 and 10.2. But after migrating to 10.3 it says Licenses are not available. Meanwhile ArcMap and ArcCatalog opens with no prob. Any idea why it may be happening?

            m_AOLicenseInitializer = new LicenseInitializer();
            if (!m_AOLicenseInitializer.InitializeApplication
                (new esriLicenseProductCode[]
                {   esriLicenseProductCode.esriLicenseProductCodeAdvanced,
                    esriLicenseProductCode.esriLicenseProductCodeBasic,
                    esriLicenseProductCode.esriLicenseProductCodeStandard,
                    esriLicenseProductCode.esriLicenseProductCodeEngine,
                    esriLicenseProductCode.esriLicenseProductCodeEngineGeoDB,
                    esriLicenseProductCode.esriLicenseProductCodeArcServer
                },new esriLicenseExtensionCode[]{}))
            {
                MessageBox.Show(m_AOLicenseInitializer.LicenseMessage());
                returnValue = false;
                m_AOLicenseInitializer.ShutdownApplication();
            }

it says that none of the licenses exits.


texts: Writing categorial rasters with color and attribute symbology info in R

texts:
I am generating landscape classification output rasters and adapted the code from WAF’s stackoverflow answer to the this question to generate GeoTIFFs with embedded standardized color and attribute information.
library(rgdal)
library(raster)

addColorTable <- function(inRstName, outRstName, rat.df){
  r<- readGDAL(inRstName)
  rat.df$color<- as.character(rat.df$color)
  rat.df$attribute<- as.character(rat.df$attribute)
  outRst <- writeGDAL(r, outRstName, type="Byte",
                      colorTable=list(rat.df$color),
                      catNames=list(rat.df$attribute), mvFlag=11L)
  return(raster(outRst))
}

# create dummy data set
r <- raster(nrow=10, ncol=10)
r[] <- 0
r[51:100] <- 1
r[3:6, 1:5] <- 2
r[1, 1] <- 3
writeRaster(r,'dummy_raster.tif',overwrite=T)

#create color and attribute table
# This defines the values, the color and the attribute
valT <- c(0,1,2,3)
colT <-  c("#FF0000", "#FF9900" ,"#99FF00","#0000FF")
attT <- c('Forest','Water body','City','Cropland')
rat.df <- data.frame(value=valT,color=colT,attribute=attT)

# apply the magic function
rnew <- addColorTable('dummy_raster.tif', 'dummy_raster_with_symbology.tif', rat.df)

plot(rnew)

The code works perfectly within R to generate the categorical land cover rasters with the proper colors (see figure above). However, no matter what I do, when I visualize the categorical raster outside of R (e.g., QGIS or ArcGIS) the raster legend includes legend elements for 0-255 raster values, even if my raster only has values between 0 and 3 (see figure below). Most importantly, element 11 in the legend is always displayed with a 0% opacity symbology despite nothing in my code specifying opacity.
Does anyone have any solutions to correctly save raster color and attribute info in R?

In the case above, the transparency of class 11 is not an issue because there are only 4 categorical classes, but in my actual classification rasters, my 11th class gets 'erased' as well.


texts: Get fields comments instead of names when viewing object by "Identify feature" instrument

texts:
I have a point layer with column NAMES like: Object_ID, Type, Link, Area an so on. Each column title has a COMMENT, like: Object identification number, Type of organization, site URL, Floor area in square meters... When I click on the point with Identify feature instrument I get a list of it's NAMES - VALUES pairs. Is there a way to replace NAMES with COMMENTS for this instrument in order to see COMMENT-VALUES pairs? May be some other instruments exists to get such kind of data for each feature while viewing the map?
Looks like I want something strange, colleagues, nevertheless... Any ideas?


texts: Giving points in the same delimited text layer multiple colours-while using point displacement

texts:
I want to give multiple colours to points within the same delimited text layer, while using point displacement, looking for advice on how best to do this.
The data are showing occurrences of a variety of species in an area. I'm displaying my points using point displacement, as there are some over-lapping co-ordinates. I would like to show a different colour for each species while still using point displacement.
I've tried splitting the data and creating a separate layer for each species, giving each species layer its own colour in that way. However, when I do this I've only been able to assign point displacement symbology to each layer individually and not to all the data at once-so I end up with overlapping points again.


texts:
As I get more into QGIS scripting, I refer to the documentation often and usually without success. I have trouble deciphering how to interpret the classes, etc.

I've been really struggling to create a New Print Layout and add things to it like the map, legend, label, etc. I have not gotten to trying to make things look nice yet, just the very basic functionality.

Let's take the QgsLayoutItemLegend class for example from the documentation:

Now, I think I understand that in order to create a legend item and add it to my Print Layout, I would do something like this...

QgsLayoutItemLegend.create(legend)

Except, I need to already have a legend item to "create" a legend in the Print Layout. But how do I create this legend item?

I don't only want to know how to create this legend item, but where in the documentation can I even find this?


texts:
Following this example and Choropleth from Mike's site I am trying to get:

image + vector + choropleth = my map

As you can see I lost the zoom capability on the vector layer. I have been trying differents things, but nothing works.

How can I fix that?


texts: How to link dataframes in ArcGIS or QGIS?

texts:
The below image shows a layout in ArcGIS 10.1 with four dataframes.  Each dataframe is of the same geographic area, although they are different images. For example, the top left image is a 1989 DOQ, the top right image is a reversed 1989 DOQ, the bottom left image is a 2012 DOQ, and the bottom right image is a reversed 2012 DOQ.

Is it possible to link the dataframes so that panning the image in
one dataframe pans the images in all of the other dataframes (i.e.
similar to how you can "Link Views" in Erdas Imagine)?
Is this possible in QGIS?


texts:
I am using QGIS 2.18.20 on Win10. I am trying to manually move my labels, and I've followed the protocol from How can I connect label and point with a line/arrow?.
I've used this protocol in another map and everything went fine, but for some reason, now whenever I try saving my changes my labels move and cluster, and end up looking like this:

What is happening?


texts: How to Sextante model builder within QGIS to convert KML to shapefile?

texts:
Are there any tutorials on how to use Sextante within QGIS? How I could use Sextante model builder in QGIS to convert KML to shapefile?

This is a routine task and Sextante would really help.

I am straggling to find standard gdal/ogr converters within the Sextante algorithms, such as shp to kml to shp to tab.


texts: Blending two overlapping DEMs using ArcMap

texts:
I have two point clouds with X,Y,Z which have been converted to DEM's.
One is more accurate than the other.
The point at which they overlap (just a 100m strip) however is very uneven.
Is there a way I can make it blend better, I have tried using focal stats but it still isn't enough.
Is there any other method I could use?
Here's a picture of the "unblended" DEM merged. I am wanting a better blend


texts: How can i create a webservice the shortest route between two POIs?

texts:
I would like to create a prototype of a webservice that can calculate the shortest route between two POIs. How can I do that ?

(I don't have any knowledge about the back-end of GIS and webmapping. So if you know any gis server or tool that can simplify this task it will be very helpful.
I prefer opensource free tools. However, A very good tool or server that worth the price will be helpful for me )


texts: Writing Python list data into CSV file in multiple columns?

texts:
I am trying to write data into csv file from python list.
I am using aws lambda for this......

My code

fields = ['SnapshotId', 'StartDate']
csvio = io.StringIO()
writer = csv.writer(csvio,delimiter=',')
writer.writerow(fields)
for x in all_snapshots:
    for snapshot in x['Snapshots']:
        snap_list=[snapshot['SnapshotId'], snapshot['StartTime'].strftime("%b %d %Y")]
        writer.writerows([snap_list])
        print(', '.join(snap_list))
s3.put_object(Body=csvio.getvalue(), ContentType='application/vnd.ms-excel', Bucket=bucket_name, Key=file_name)

This code is pushing csv file to s3 but all the list values are coming into single csv column with comma seperation. I want to put this into multiple columns.

THis is the csv file generated:

SnapshotId,StartDate

snap-aaaaaaaa,May 25 2016

snap-aaaaaaaaa,Jul 14 2016

This is into single column, but I want seperate columns for snappshot id and start date.


texts: Geotagging photos with Point Shape via Time Stamp

texts:
I have waypoints which I marked with a GPS device - unfortunately no track/route. And I have photos which I took at the same day.
Now I want to interpolate the points to line and use this line to geotag my photos.

With QGIS and the plugin "point to path" it's easy to get a line shape. This line shape can be saved as a gpx-track.
With the Freeware "GeoSetter" it's easy to geotag photos with a gpx-track.

But when I try this (Point to Path -> save line shape as gpx -> geotag photos with gpx), GeoSetter complains that the gpx-file doesn't have coordinates with time stamps.

So, how do I get the time stamp of the waypoints into the gpx-track?


texts: How to get rid of unsharp picture markers?

texts:
I have made a symbol in Adobe Illustrator to use for a custom picture marker. The thing is when I select these markers and put them in my map the get blurry. I have tried several image formats such as .png .jpg .emf and .bmp. But no succes.

I was looking if there is a way to use a .svg file, but with no result.


texts:
I am new to QGIS WMS and trying to add imagery from Six Maps.

https://maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_Imagery/MapServer

I have tried the steps in here also adding as an ArcGisFeatureServer -> New connection QGIS: Connect to a REST service however this didn't work either.

I get the error message

Failed to retrieve service capabilities:
Parsing error: illegal value


texts: Calculating Optimal Geohash Precision from Bounding Box

texts:
I'm using Elasticsearch's GeoHash grid Aggregation to plot clusters on a map (using Leaflet). I understand that for larger areas a lower precision setting should be used to limit the number of buckets created/returned.

How should I determine the appropriate precision value to request?

Is there a standard or recommended formula for calculating the optimal precision based on a bounding box and/or zoom level? Is it better to just map zoom levels to precision values (I know that's probably easiest.)


texts: Apply different symbology layers to different dataframes with arcpy mapping?

texts:
I have one shapefile, which I want to create a map of with 3 different data frames (showing three different attributes of the shapefile) on one map. I have 3 separate symbology layers files saved to be applied to each dataframe. I have got as far as loading in the shapefile to each data frame, and can apply the same symbology layer to each, but can't get different symbology layers to go to different dataframes. I want this to be automated as eventually I want to loop this over a large number of shapefiles to make multiple maps, but I just can't work out this beginning step.

So far I have:

import arcpy,

shape = arcpy.mapping.Layer("myshapefile.shp")
symbologyLayer =("mysymbology.lyr")

for df in arcpy.mapping.ListDataFrames(mxd):
 arcpy.mapping.AddLayer(df, shape, "AUTO_ARRANGE")
 arcpy.ApplySymbologyFromLayer_management (shape, symbologyLayer)


texts: OpenLayers - How to display selected features as graphic markers?

texts:
I'm attempting to render a point feature as a graphic.  When that feature is selected I'd like to display a green indicator above the graphic.  It seems backgroundGraphic and externalGraphic are perfect for this...but I can't get it right.

    function getStyles()
    {
        var unselectedMarkerOptions = {
            pointRadius : 24,
            strokeWidth : 1,
            fillOpacity : 1,
            graphicHeight : 48,
            graphicWidth : 48,
            externalGraphic : "http://usgvdcalix2:83/openlayers/img/blank.png",
            backgroundWidth : 48,
            backgroundHeight : 48,
            backgroundGraphic : "http://usgvdcalix2:83/openlayers/img/marker-blue.png"
        };

        var selectedMarkerOptions = {
            pointRadius : 24,
            strokeWidth : 1,
            fillOpacity : 1,
            graphicHeight : 8,
            graphicWidth : 8,
            externalGraphic : "http://usgvdcalix2:83/openlayers/img/greendot.png",
            backgroundWidth : 48,
            backgroundHeight : 48,
            backgroundGraphic : "http://usgvdcalix2:83/openlayers/img/marker.png"
        };

        return {
            selected : new ol.Style(ol.Util.applyDefaults(selectedMarkerOptions, ol.Feature.Vector.style["default"])),
            unselected : new ol.Style(ol.Util.applyDefaults(unselectedMarkerOptions, ol.Feature.Vector.style["default"]))
        };
    }

I add these styles to a stylemap and associate it with a vector layer but it doesn't behave well:

var styles = getStyles ();
var styleMap = new ol.StyleMap({
    "default": styles.unselected,
    "select": styles.selected
});

var featureLayer = new ol.Layer.Vector ( "Features", {
    styleMap : styleMap,
    renderers : ol.Layer.Vector.prototype.renderers
} );

It requires a transparent image (blank.png).
User must click the externalGraphic meaning the 8x8 "greendot.png".
Setting backgroundGraphicZIndex to a large value to make it a foreground graphic didn't work.

The OL select-feature-multilayer example works well.  But I'm looking for a tiny green indicator to show that a feature is selected with the ability to click the feature and toggle the indicator.  How is this possible?


texts: Display the resulting value of the sql.execute query into a

texts:
I need to display the total number of rows(sum) in a DIV element. The number would express the sum of the data being toggled using the following function:

var LayerActions = {
  extortion: function(){

    sublayers[0].set({
        sql: "SELECT * FROM table1 where type <> 0",
    });
    return true;
  },
}


texts: Replacing data source of ArcMap Service Layer using ArcPy?

texts:
I have many map documents that point to a hosted imagery service layer that has changed names.

I would like to write an ArcPy script that will replace the data sources of service layers that have changed names.  So far I have tried the following from the Python window, without success (assuming the only layer in the MXD is the service layer):

mxd = arcpy.mapping.MapDocument("CURRENT")
lyr = arcpy.mapping.ListLayers(mxd)[0]
sourceURL = lyr.serviceProperties["URL"]
arcpy.AddMessage("Source URL: {0}".format(sourceURL))
destURL = "https:\path\to\image\service"
arcpy.AddMessage("Destination URL: {0}".format(destURL))

lyr.findAndReplaceWorkspacePath(sourceURL, destURL)

arcpy.RefreshTOC()
arcpy.RefreshActiveView()
mxd.save()


texts:
I am trying to import a .shp file and CSV file so that I can join data. I have written below code in Python console within QGIS.  For Layer.isvalid(), I am getting False. I am using QGIS 3 and on Windows 7. I don't know where I am doing mistake. Since I am using Python console inside QGIS should I need to set enviorment variables?

layer = QgsVectorLayer("file:///C:Test/t1_2013_06_tract/t1_2013_06_tract.shp","BuildingLayer","ogr")

layer.isValid()


texts: Shapely LineString contains method not showing contained points

texts:
I have some GeoJSON files that I converted from some KML someone gave me to analyze.

I have used GPSBabel to convert them from KML to GeoJSON. The reason for doing so is that I am able to easily process them using python 3.5 libraries like json and shapely (my environment does not work with packages like fiona for version incompatibilities).

A sample command is: gpsbabel -w -t -i kml -f input.kml -o geojson -F output.json

The KML contains various <Polygon> each one with one <LinearRing> inside.
The GeoJSON contains various geometry coordinates of type called LineString (this comes from the GPSBabel automatic conversion).
Each point in the KML/GeoJSON is a 3 fields array which should be [longitude,latitude,altitude] (according to this: https://developers.google.com/kml/documentation/kmlreference#linearring).

I've written the following function to establish if an input Point (Shapely geometry) is contained in the LineStrings (Shapely geometry) found in the input:

from shapely.geometry import Point, LineString

def poly_contains(point, poly_js):
    for feature in poly_js['features']:
        lineString = LineString(feature['geometry']['coordinates'])
        if lineString.contains(point):
            return True
    return False

Each input point is created from a CSV containing the coordinate pairs with something like this: Point(row['LatitudeId'], row['LongitudeId'])

When applying the previous function to all these points I can not find any point contained in those shapes, but I know for sure some of those points are indeed contained!

I thought the reason could be the first and the last point of a LineString don't match (i.e. the "polyline" is not a "polygon"), but the following piece of code makes sure I am dealing with "polygons" defined with matching first and last element:

def check_file(file_content):
    for feature in file_content['features']:
        lineString = LineString(feature['geometry']['coordinates'])
        if lineString.coords[-1] != lineString.coords[0]:
            print("Issue!")

When applying the previous function to these GeoJSONs like this:

import json
with open('.../.../myfile.json') as f:
    stuff = json.load(f)

check_file(stuff)

then I don't see any output like Issue!, so all these LineStrings should represent "polygons" right?

I don't have many details about those projections and stuff like that, but I am assuming this is all Web Mercator in WGS84 coordinate system as this data is used to be represented on Google Maps (both the original KMLs and the points from the CSV that I need to check).

What's wrong with what I am doing? Why that function contains applied to LinearString and Point does not work as expected?

Edit:

I may be confusing WGS84 and Web Mercator (EPSG:3857), I am not sure. However I have tried to project all those shapes (the points from the CSV and the LineString from the KMLs/GeoJSONs) using pyproj, then wrapping those shapes in Polygon instead of LineString like this:

from shapely.geometry import Point, Polygon
from pyproj import Proj
p_web = Proj(init='EPSG:3857')

# applying 2D Web Mercator projection
def poly_contains_v2(point, poly_js):
    point_proj = Point(p_web(point.coords[0][0], point.coords[0][1]))
    for feature in poly_js['features']:
        new_coords = []
        for ring in feature['geometry']['coordinates']:
            x2, y2 = p_web(ring[0], ring[1])
            new_coords.append([x2, y2])
        polygon = Polygon(new_coords)
        if polygon.contains(point_proj):
            return True
    return False

but still I can not find any point included in those "polygons". I have also tried to apply just a Polygon wrap instead of a LineString (with no projection) with no luck.


texts: GPSbabel: how do I define the layout of my input data

texts:
I have collected GPS data on a journey by extracting it every second from the $GPRMC NMEA sentence, with each reading saved in the form:

time,latitude,longitude,speed,direction

a typical couple of lines are.

163044,5050.8427,-00049.8211,047.4,107.6
163045,5050.8387,-00049.8012,047.2,106.5

Note: I am using Linux, and a Processing sketch to extract the data and save it in a file made unique by including the date and time of its creation in its name, e.g. posns20150321131751.txt

I intend to plot this as a trace on an OpenStreetMap which requires data in osm format.

This is my Linux command line and error message:

[Harry@localhost]~/Gpsbabel% gpsbabel -i xcsv,style=/home/Harry/Gpsbabel/myidata.style -f /GPStalk/posns20150321131751.txt -o osm ~/GPStalk/posn01.osm
CSV_UTIL : xcsv style "%H%M%S" is missing default.

My style file myidata.style is:

FIELD_DELIMITER    COMMA
IFIELD GMT_TIME,"","%H%M%S"
IFIELD LAT_NMEA
IFIELD LON_NMEA
IFIELD IGNORE
IFIELD IGNORE

I find the documentation in the GPSBabel documentation difficult to follow, so my question is, please, where have I gone wrong? In particular, what is wrong with my style file myidata.style?


texts: How to define layer's URL programmatically? (ESRI Silverlight API)

texts:
I have defined a simple tiled map service layer:

<esri:ArcGISTiledMapServiceLayer x:Name="BaseMap"
ID="BaseMap" Opacity="50" Initialized="MapServiceLayer_Initialized"
InitializationFailed="Layer_InitializationFailed" />

Then I assign its Url property programmatically:

BaseMap.Url = App._componentDict["BaseMapRestUrl"];

The app bombs with this message, and the entire Silverlight app goes blank:

Layer 'BaseMap' is currently unavailable. Error" Value cannot be null. Parameter name: Url

Is there a particular way to go about this that I'm missing? I suspect I'm setting the Url too late in the lifecycle.


texts: At the request of my superiors, I'm having to make 1, sometimes multiple changes, to Web App Builder Developer sites throughout the day that include (searches, pop-ups, custom widgets, etc).  Upon reloading the new site files to IIS, users sometimes still see the previous settings in the sites.  I have to get them to actually clear their browser cache in order to see the new changes.  Is there any way around this??  I'm hoping there is some setting in IIS that may help with this..


texts:
The output I get from QGIS 7.0.1 for an ASCII DEM file doesn't seem to be a correct format.  It differs from the format I get when I download a DEM file from webGIS ( http://www.webgis.com/terr_us75m.html) .

I have an application that will accept an ASCII DEM file as input.  When I input the DEM file created by QGIS it says the format is invalid.  When I input the DEM file I downloaded from webGIS is fine.

The start of the QGIS DEM file looks like this:

         gmted_mea075_Contour-02.dem                                                                     -1223458.8333  4712 3.7222

The start of the downloaded DEM file looks like this:

SPANAWAY, WA-24000  LAT::   47  0  0.0000 N LONG:: -122 22 30.0000 W SCALE:: 24000   (C) www.webmet.com

Has there been a change in the ASCII DEM format?

My input application is from 2003.

I create the DEM file:  Raster > Analysis > DEM (Terrain models)
For the file type I select:  USGS Optional ASCII DEM(*.dem *,DEM)


texts: How do I transfer specific attributes from a point feature class to a polyline feature class

texts: I have a point feature class that I have used to create stream segments of interest. Now I'm trying to transfer some of the attributes from the points feature class to the new polyline feature class. Transferring all of the attributes would be ok, but I'm not interested in all of them. I've tried the Attribute Transfer Mapping function under the Spatial Adjustment Toolbar and think that's what I need but have evidently not been feeding the tool the right info. I'm using ArcMap 10.0. Thanks!


texts: What CRS to use for UTM 36K Georeferencing?

texts:
I’m trying to georeference some old scanned maps in QGIS, but I’m not sure what CRS I should be using.

The legends on the maps say:

  Numbered lines indicate the 1,000 metre Universal Transverse Mercator
  grid, zone 36K

  Modified Clarke 1880 spheroid

WGS 84 / UTM zone 36S (EPSG:32736) seems to be giving me the best results, but it’s still not satisfactory. See particularly the tile overlaps and streams in the image below. If I choose WGS 83 (EPSG:4326) my georeferenced images disappear from view.

I’ve been using “Linear” for the transformation type and the resampling method in the georeferencer. Might it help to use other settings here?


texts:
I'm  using ol3 in a  mobile  hybrid application using sencha touch2  , so I start by drawing a polygon inside a map , after validating this drawing I  redirect the user to an other form in order to fill the alphanumeric  informations of that polygon , My problem is when I want to get back  to  my map  the layers disappear and the only way to make them get back is to change the device Orientation.  so I'm looking for code snippet that can reload (refresh) my layers or the whole map

I tried
    map.dispatchChangeEvent();

or vectorSource.dispatchChangeEvent();

but nothing works for me

I find the solution is to update the map size  like  this : map.updateSize();

I


texts: Infering course given origin and destination coordinates

texts:
I have two points. A vehicle is at point A and I have the coordinates at point A. The vehicle moves to point B and I know the coordinates as well.

I would like to find the angle of movement i.e., 0,180,90,45,360 degrees etc. With the angle, I can tell in what direction its headed relative to origin.

Given that I know the coordinates at A and B, how would I calculate the degree between point A and B.


texts: Split polyline where it intersects with other lines in MapInfo

texts:
I have a polyline (bus route) which I want to split where it intersects with other lines (not polylines) I have created.

I'm happy to initially do it manually to see how it works, although I would ideally prefer to do it all entirely in MapBasic due to the amount of data I have to get through.


texts: Draw Bounds with coordinates in QGIS from CSV

texts:
I have a CSV file that contains 4 coordinates on each line, 1 min latitude, 1 max latitude, 1 min longitude and 1 max longitude. These 4 points delimit a particular region(they are the bounds) I would like to draw for each row of my CSV, the region corresponding to the 4 points. How can I do that?

To be more clear this is what i want(done with OpenStreetMap):

This is the csv:

In particular NORTH,SOUTH,WEST and EAST are the 4 coordinates that define my bounds


texts: How to set the fill extrusion height in Mapbox Studio

texts:
in Mapbox Studio, I'm trying to set the height of an extrusion layer to a value that is mapped to the identity values of the layer property but I cannot find the option.

My layer data have the property upper_val_m that contains the height
in meters of each polygon.

In Mapbox Studio I have the option to map the height value to the data value "exponential", "interval" or "categorical" but none of those seems to allow an "identity" mapping directly to the data value.

What did I miss?


texts: Change printer for a specific map document from Arcpy

texts: Is it possible to change the printer for a particular map document from ArcPy? The mxd is saved with the default printer settings, and I know with the PrintMap() function one can specify the printer to use. In my case, I do not want to print the map, just change the printer.


texts: Rectangle polygon is displayed as triangle in QGIS?

texts:
I am wondering why a rectangle polygon that i have created on an official german surveying website (downloaded as compressed shape in a zip-file) is displayed as a triangle when loaded in qgis. mighty this have somethting to do with the chosen refence system (wgs84)?

at first this problem could temporarely be solved by starting the editing mode of this layer, but when ending editing the displayed rectangle is shown as a triangle again.


texts: Displaying rotated grid data in Leaflet with GeoServer WMS

texts:
I am trying to display rotated grid data from NetCDF in Leaflet.
When I try to use the "WMS angle parameter" (https://docs.geoserver.org/latest/en/user/services/wms/vendor.html#angle)
I got some strange results like this:

Any idea why this happens?
I already have seen that there is the possibility to create a custom projection which allows to add rotation parameters, but I could not figure out how the correct usage is with using my own projection. Currently I am using (EPSG:25832)
(Defining coordinate reference system with rotation in GeoServer?)

Would be nice to find a solution to display the rotation without manipulating the raw NetCDF data.


texts:
I used QGIS for a very long time.
In the last version, there seems be a new sound (buzz) when the program is run.
Does anybody know how we can remove this sound?


texts: Copying polygon geometry in ArcPy using old style cursor?

texts:
What I need to do is to create a copy of a polygon shapefile to a new shapefile with slightly different attributes.

Where I got stuck is in copying the geometry - Geometry doesn't have to change at all in the new shapefile. The actual 'copying' of the polygon geometry is not the problem.
I loop through the input shapefile:

> for row in rows:
>     partnum = 0
>     feat=row.getValue("SHAPE")

And then I enter this shape in the output feature class:

outFeat.Shape=feat

Problem is, I also want to calculate the Area, Perimeter and centroid of the polygon and I can't seem to 'cast' the feat variable to a polygon so I can then get the polygon properties (area, centroid, etc). Can you do that somehow?

One of the things I tried was to create the polygon from scratch, i.e. loop though the vertices of the input polygon and create a new polygon/array to use for the new polygon. This method apart from being IMO over the top, it also had the problems with the multi-part polygons. My polygons ended up corrupted - although I followed the steps outlined here

I didn't dwell too much with it since what I did instead is create a dummy Polygon which I then passed to it the current geometry value in the cursor.

# Create a dummy polygon
            dummyArray=arcpy.Array()
            dummyPoint = arcpy.Point()
            coordList = [[[1,2], [2,4], [3,7]],[[6,8], [5,7], [7,2], [9,5]]]
            for feature in coordList:
                for coordPair in feature:
                    dummyPoint.X = coordPair[0]
                    dummyPoint.Y = coordPair[1]
                    dummyArray.add(dummyPoint)
                dummyArray.add(dummyArray.getObject(0))

            parcel=arcpy.Polygon(dummyArray)
            parcel=feat ### --> This is where I replace the dummy geometry value with the current value in the loop
            parcel_c=arcpy.Point()
            parcel_c=parcel.centroid # I can now get to the polygon properties
            c_x=parcel_c.X
            c_y=parcel_c.Y

Although this works fine, I am not really happy with it as its a bit of a hacking and -as mentioned in the beginning- I think I am missing something vital!

Any insights?


texts:
I'd like to generate a report about MXD layers for a handful of MXDs. I want the report to display their name, workspace path and scale range settings.

I can get the first two with ArcPy but I didn't find a way to get the scale ranges. Can this be done with ArcPy? Or perhaps there's some other way to generate the type of report I'm looking for?

I found this but it doesn't give me scale range info: http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=A910AB18-1422-2418-3418-3885D388EF60

Here's the code I've started with:

import os, arcpy

directory = r"D:\MapTemplates\Map_Templates"
for root, dirs, files in os.walk(directory):
    for myFile in files:
            fileExt = os.path.splitext(myFile)[1]
            if (fileExt == ".mxd"):
                fullPath = os.path.join(root, myFile)
                print myFile

                myMap = arcpy.mapping.MapDocument(fullPath)
                layers = arcpy.mapping.ListLayers(myMap)
                for layer in layers:
                    wsPath = ""
                    try:
                        wsPath = layer.workspacePath
                    except:
                        pass
                    if wsPath:
                        print "    " + layer.name + " : " + wsPath


texts:
I am at the end of my rope here. This cursor runs fine for one iteration and then stops. It doesn't throw an error, it just completes. There are 3,000 rows in my table, so I don't know why it would do this. Does the Generate Near Table function somehow stop the cursor?

import arcpy
rows = arcpy.SearchCursor("mergedVMS")
trip_id = 0
for row in rows:
    if row.trip_id == trip_id:
        continue
    else:
        trip_id = int(row.trip_id)
        arcpy.SelectLayerByAttribute_management("mergedTrawls", "NEW_SELECTION", '"trip_id" = ' + str(trip_id))
        arcpy.SelectLayerByAttribute_management("mergedVMS", "NEW_SELECTION", '"trip_id" = ' + str(trip_id))
        arcpy.GenerateNearTable_analysis("mergedVMS", "mergedTrawls", "C:\\Desktop\\Feb25\\Project.gdb\\pointdist" + str(trip_id))
del row, rows

Per dmahr's answer, my new script that works is as follows:

import arcpy
rows = arcpy.SearchCursor("mergedVMS2")
trip_id = 0
for row in rows:
    if row.trip_id == trip_id:
        continue
    else:
        trip_id = int(row.trip_id)
        arcpy.SelectLayerByAttribute_management("mergedTrawls", "NEW_SELECTION", '"trip_id" = ' + str(trip_id))
        arcpy.SelectLayerByAttribute_management("mergedVMS", "NEW_SELECTION", '"trip_id" = ' + str(trip_id))
        arcpy.GenerateNearTable_analysis("mergedVMS", "mergedTrawls", "C:\\Desktop\\Feb25\\Project.gdb\\pointdist" + str(trip_id))
        arcpy.SelectLayerByAttribute_management("mergedTrawls", "CLEAR_SELECTION")
        arcpy.SelectLayerByAttribute_management("mergedVMS", "CLEAR_SELECTION")
del row, rows

I have simply duplicated my mergedVMS layer and run the cursor on that (trip_id is my unique value). I also added two clear selection lines. It is currently running exactly as expected.


texts: gpTask: ERROR 000735 Value is required (Input KML, DataFile)

texts:
I have a Python script tool running perfectly on the ArcGIS server. Trying to run the tool in a web application (JS API 4.6) with the same input KML-file (published on a publicly available webserver). I get always the error ERROR 000735 Value is required.

The input KML is defined as a DataFile by using the URL. In the console log, I can trace the DataFile object - the URL is correct and the object defined. The tool requires a GPDataFile as input.

var toolUrl="https://...";
function runTool() {
    var inputFile = new DataFile({url:"http://.../test.kml",});
    var params = {"inputFile":inputFile};
    var gpTask = new Geoprocessor(toolUrl);
    gpTask.execute(params).then(drawResultData);
  }


texts:
I want to import a DWG file into QGIS (3.8.1). I do this via Project > Import/Export > Import Layers from DWG/DXF.

Everything works fine, the imported data looks almost like the original DWG - with one difference. All CAD layers which only contain points are not imported. I first thought, this might be because I ticked the option to "merge layers" (so i got only 4 layers "hatches", "texts", "lines" and "polylines" as result).
Importing the data to another geopackage without merging layers did not change the issue - in the layer overview in QGIS of the DWG the point layer exists, and i have checked it to import it's data - but it does not show up in QGIS in the resulting group.
Any ideas why this is happening and how to solve this?


texts: Custom locator not matching to entry with best score using ArcGIS Pro

texts:
I created a custom locator using the Create Locator tool to create individual county geocoders, then used the Create Composite Address Locator to combine them. When I geocoded a CSV file,  I had about 19% that were "Tied".

As I started using the Rematch Address tool to fix them, I noticed that there was actually a 100% match for most of them, but Arc defaulted to the one with the worst score.
Example:

Of those that were "Tied", Arc automatically displayed the points with a 60% match even though there is a 100% match.
When I tried looking for a solution in help, all I found was that it should've matched to the best score.

https://pro.arcgis.com/en/pro-app/latest/help/data/geocoding/tips-for-improving-geocoding-quality.htm
Seeing how there are over 86,000 instances of this, I'd prefer not to go through every entry to correct. Is there a way to just rematch all of them to the best score?


texts:
I have a PNG map with just black and white colors, where black is the one pixel border. I want this in vector, or want an easy way to draw 90 degree polygons.

I tried some things from forums but I don't have a solution yet.


texts: Update a count of the number of polygons which intesect with polygons in the update table

texts:
I have a table with 1km squares as polygons and I want to update a column with the number of spatial intersects in another table of polygons. Something like this [if count could be used here];

Update KmStats
SET RecCount = COUNT(*)
From KmStats
INNER JOIN SP_Records
   On KmStats.SP_geometry.STIntersects(SP_Records.SP_GEOMETRY) = 1

All the polygons in SP_Records are aligned to the same grid and are either small enough to fit in a 1km sq or are 1km or 2km sq polygons.


texts: User locking postgres/postigs table without privileges

texts:
I have a PostgreSQL 8.4.4 database with PostGIS 1.4.2 and I am experiencing a problem where an external GIS system (using embedded CadCorp SIS 7.0) is connecting to the database and taking Access Share Locks on a large number of tables which, as far as I can understand, it should not have privileges to access. The affected tables are all sources for views that are registered in geometry_columns, across a number of schemas. All other tables, whether registered in geometry_columns or not, with or without a geometry column, are unaffected (except those explicitly required by the map project the external system is using). This is frustrating because the connections stay open for a large period of time and prevent a lot of maintenance work that I might like to do.

For example I have a schema:

CREATE SCHEMA my_schema AUTHORIZATION user_one;
GRANT ALL ON SCHEMA my_schema TO user_one;

And within that a table:

CREATE TABLE my_schema.my_table...;
ALTER TABLE my_schema.my_table OWNER TO user_one;
GRANT ALL ON TABLE my_schema.my_table TO user_one;

And a view:

CREATE OR REPLACE VIEW my_schema.vw_my_view AS SELECT... FROM my_schema.my_table;
ALTER TABLE my_schema.vw_my_view OWNER TO superuser;
GRANT ALL ON TABLE my_schema.vw_my_view TO superuser;
GRANT SELECT ON TABLE my_schema.vw_my_view TO user_one;

The roles user_one and user_two are independent but I find that user_two is able to put an Access Share Lock on my_schema.my_table. If I connect to the database as user_two I am denied access to my_schema.my_table as I would expect.

Is there something I'm not doing with the privileges or roles that would prevent this? Why is it only those tables that are the source for views that are registered in geometry_columns that are affected?


texts: ArcGIS Identity tool giving ERROR 999999

texts:
When I run the Identity tool in ArcGIS software I get this error.

How can I fix it?


texts: ArcPy - error with updateRow() using UpdateCursor

texts:
I am trying to use the UpdateCursor function to run through each field in each row and to compile a string which is then added to the AbvField.
I'm trying to do this in the python window, rather than as a standalone script.

I get the error:
Runtime error : Row: Field updateRow does not exist

I have included all the code except for the subs dict, which is simply a dictionary of species names (keys) and abbreviations (values) with which the script creates the string for the AbvField.
Everything works fine except for the updating the field.

TheShapefile="Primary_biodiversity_layer"
AbvField = "Sp1"

Rows=arcpy.UpdateCursor(TheShapefile)

Fields = arcpy.ListFields(TheShapefile)

SpList = []
newDict = {}

for TheRow in Rows:
        RString = ""
    for TheField in Fields:
        TheFieldName=TheField.name
        TheValue=TheRow.getValue(TheFieldName)
                if TheValue == "":
                    pass
                elif TheValue not in subs.keys():
                    pass
                else:
                    RString += subs[TheValue] + ", "
                    if TheValue not in newDict.keys():
                        newDict[TheValue] = subs[TheValue]
                    else:
                        pass
        TheRow.setValue(AbvField, RString)
        TheRow.updateRow(TheRow)
        del TheRow


texts: Completing dashed line shapefile feature that was created with ArcScan?

texts:
Using ArcScan to vectorize a raster floodplain boundary map, I have selected the area, run ArcScan in the area, to output a Shapefile of the line. ArcScan worked to create the line, and it is a dashed line, so it created a dashed line. However, I need to connect the line. I have used the editor, continue line and connected them at the end points, but I was wondering if there is any way I can automate the joining of the dashed line.

I have tried:

Advanced Editing, Extend Tool, but it only works if the lines are perpendicular to one another.

Editor, Buffer, did not work

I am thinking if I could save the ShapeFile, then re-rasterize it, use ArcScan again on the cleaner dashed line, and then turn on the ArcScan/Vectorization Settings, Gap Closure Tolerance and then it would close the gaps.

I will post images below.


texts:
I use foursquare because it contains a lot of POIs.
Which database do you use when you need a point of interest in North America?


texts: Cannot save any geoprocessing result with ArcPy?

texts:
I have a problem to get the code below running as a python script tool in ArcGis 10.2.2. I always get this error Message:
AttributeError: ResultObject: Get attribute save not found
The tool help clearly says that there is a .save function. What am I doing wrong?

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

inRaster = arcpy.GetParameterAsText(0)
DGM = arcpy.Raster(inRaster)

# Process: Fill
Fill_DGM = arcpy.gp.Fill_sa(DGM, "Fill_DGM", "")
Fill_DGM.save("F:/Desktop/test/Filled")


texts:
I know this question was asked several times and I searched plenty of past answers and tried plenty of variations but I couldn't transform my coordinates correctly.
My coordinates (stored in a .csv) are given in WGS84 Lon/Lat decimal degree format.
For example I have the coordinate Lat: 38.906 Lon: -120.71
This should output the UTM coordinates: 698568.503, 4308838.330
But if I do (following this instructions):
cord.dec = SpatialPoints(cbind(mycsv$Longitude,
                               - mycsv$Latitude),
                         proj4string=CRS("+proj=longlat"))
cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32610"))

# add coordinates to UTM list
utm$X=cord.UTM@coords[,1]
utm$Y=cord.UTM@coords[,2]

I get: 698714, -4308793 which is obviously not correct.
I tried another instruction on this topic. But if I do:
lonlat.coordinates <- data.frame(ID = 1:length(mycsv$Latitude),
                      Y_LAT = mycsv$Latitude,
                      X_LON = mycsv$Longitude)

coordinates(lonlat.coo) <- c('X_LON', 'Y_LAT')

I get:

Error in [.data.frame(object, , value) : undefined columns selected

Can anyone please help me?
Best wishes


texts: Intelligence Tools don't show up in ArcGIS Pro

texts: I need to use the "Points to Track" tool in ArcGIS Pro, but it does not show up when I search for it in the geoprocessing toolbar. I searched online and found that it is in Intelligence Tools. For whatever reason, the Intelligence Toolbox isn't available in my version of ArcGIS Pro. Do any of you know what the problem is, and howI can resolve it?


texts: Running QGIS plugin developed for old version in newer version?

texts:
My company developed a non-public plugin to QGIS 2.6 in 2013. Now we have to run this plugin again.

Is it possible to make it run in QGIS 2.18?


texts: What is the relative popularity of ESRI GDB files vs GDB over RDBMS for custom web apps?

texts:
Our organization is interested in knowing the frequency with which people who are using ESRI geodatabases choose GDB files over GDB in relational databases for serving web applications.  Opinions and anecdotal evidence are welcome, as well as any hard data on this topic that anyone might have.

Thanks!


texts: How to meassure the length of common edges between polygons?

texts:
I have a shapefile that contains several polygons which in fact represent the limits of several counties. My task is to find each length of each shared border. So in other words I want to know how many km of border share the A county with the B county and so on. Of course I am looking for a tool to do this automatically. Measuring them manually would take ages. I envision some sort of matrix, but not sure how to do that in ArcGIS.

I use ArcGIS Desktop 10.0 (ArcEditor) .

10x


texts: Issues with Selection and Viewing of Attribute Table using QGIS and SQL Server Express 2008 R2

texts: I am having issues with Selection using QGIS 2.0.1. Dufour. When I select one point, all of the points get highlighted. Attribute Table dialog shows Errors for all of the attributes. It indicates that Features total:12691, filtered:12691, selected: 1. Has anyone else experienced this issue?


texts: Why is QGIS only one label rendering at a time?

texts:
I'm new to QGIS, but I'm experiencing what I think is a bug.

When I load the qgis_sample_data landcover raster and lakes vector layers, then enable labels for NAMES on that layer, I can only ever view one label at a time.

I've disabled simplified geometry and enabled "show all labels for this layer (including colliding labels) still to no avail.

I've attached an imgur album of screenshots to demonstrate the problem - the only thing I'm changing in those screenshots is the zoom level or pan.

Using QGIS 2.8.3 for Mac OS X Yosemite.


texts: How to merge polygons using GeoDataFrame?

texts:
I would like to merge polygons (parcels) which are inside the white road (which is a hole) -- see the picture.
The way I'm doing this is certainly a bad way because later, in my script, I have this error:
  AttributeError: 'GeoSeries' object has no attribute '_geom'

pa_fu = parcelles.unary_union
voirie = commune.difference(pa_fu) # commune is a big polygon which contains all the geometries

public = pd.DataFrame()
public['geometry'] = voirie
public['statut'] = ['public'] * public.size

prive = pd.DataFrame()
prive['geometry'] = pa_fu
prive['statut'] = ['prive'] * prive.size

decoupage = GeoDataFrame(pd.concat([public, prive]))
decoupage.crs=commune.crs # on fixe le système de projection identique à commune
cases_decoupage = overlay(cases, decoupage, how="intersection")

How can I merge these parcels and create a GeoDataFrame containing merged parcels and the white hole as a polygon?


texts: Is there a service or db to download or convert Postal Codes to Geo-Shapes / Polygons?

texts:
I'm looking for a service that converts Postal Code of Major European Cities to Geo Shapes (List of Lat/Lng) pairs like when you go to Google Maps and type for examples 10557 Berlin, and it draws boundaries around the are for that Geo-Code.

I found some Hacky ways to use the official Maps for this, But I don't wanna rely on that. is there a paid service or downloadable files that has such data.

Given that I don't want Postal Code to (1 Center Point Lat/Lng) but instead to a lit of points to form the shape for that Postal Code.


texts: Boundless OpenGeo Explorer fails with Python errors

texts:
I am running this stack on a Windows 2012 R2 Datacenter.

PostgreSQL 9.3.10
PostGIS 2.2.0
PGAdmin III 1.18.1
Geoserver 2.7.3
QGIS 2.12.0
Boundless OpenGeo Explorer 0.6.2

Everything runs rather fine, except Creation of Geoserver Catalogs and publishing. Catalogs can't be created after the Default one, but the publishing works after a fashion, i e table names can't be changed for instance. The following two log excerpts show the Catalog error and the publishing error.

2015-11-12T20:27:28 1 Traceback (most recent call last):
 File "C:/Users/Administrator/.qgis2/python/plugins\opengeo\gui\gsexploreritems.py", line 231, in <lambda>
 createCatalogAction.triggered.connect(lambda: self.addGeoServerCatalog(explorer))
 File "C:/Users/Administrator/.qgis2/python/plugins\opengeo\gui\gsexploreritems.py", line 235, in addGeoServerCatalog
 dlg = DefineCatalogDialog(explorer)
 File "C:/Users/Administrator/.qgis2/python/plugins\opengeo\gui\dialogs\catalogdialog.py", line 14, in __init__
 self.initGui()
 File "C:/Users/Administrator/.qgis2/python/plugins\opengeo\gui\dialogs\catalogdialog.py", line 92, in initGui
 self.certWidget = QgsAuthConfigSelect( keypasssupported = False)
 AttributeError: 'keypasssupported' is not a Qt property or a signal

2015-11-12T20:28:36 1 Traceback (most recent call last):
 File "C:/Users/Administrator/.qgis2/python/plugins\opengeo\gui\qgsexploreritems.py", line 59, in <lambda>
 publishProjectAction.triggered.connect(lambda: self.publishProject(tree, explorer))
 File "C:/Users/Administrator/.qgis2/python/plugins\opengeo\gui\qgsexploreritems.py", line 102, in publishProject
 tree.findAllItems(catalog)[0].refreshContent(explorer)
 AttributeError: 'NoneType' object has no attribute 'refreshContent'


texts: CGIAR SRTM to xyz using GRASS GIS, with sampling

texts:
I have used CGIAR interface to select the tiles, like so:

This resulted in a bunch of *.zip-files, and I have extracted them all into
a directory.
Next, I followed the instructions as https://grasswiki.osgeo.org/wiki/HOWTO_import_SRTM_elevation_data to turn everything into mosaic,
like so:

 #create mosaik (optionally reproject on the fly with -t_srs)
 gdalwarp *.tif srtm_mosaik.tif
 # import
 r.in.gdal input=srtm_mosaik output=srtm_mosaik location=srtm_mosaic_cgiar

Finally, I "connect" to the location srtm_mosaic_cgiar using grass-gui and export the tiles as xyz elevation data, using r.out.xyz (I use the File->Export... drop-down menu).
Now my question is, since this operation appears to be hard on my laptop,
and I can imagine the elevation data of Europe (as I have selected)
is going to be huge, how to "restrict" r.out.xyz to sample (somehow) the
exported data, or how to ascertain beforehand the size of the output?

My final goal is rather simple. I would like to have a file with x,y,z, where x,y refer to a "locus" somehwere in Europe, all the x,ys are mutually sufficiently apart (say, 300-500km), and z means the "elevation"
of the respective point x,y. If there is a misconception on my part at this
stage already, please feel free to correct me. I am new to GIS.


texts: Color-correcting/matching 4-band PlanetScope imagery

texts:
I've been trying to color correct/match 4-band (R, G, B, NIR), 3m resolution Planetscope surface reflectance imagery in order to create a mosaic for land-cover classification. At this point, I believe I've tried all possible mosaic and raster catalog color matching and color balancing options in both ArcMap 10.8 and ArcGIS Pro.
Are there some preprocessing steps that need to happen first, or is it possible there's too much variation in the imagery to create a workable mosaic?


texts: Merging isolines into polygons to produce isochrone map in QGIS

texts:
There is a set of isolines produced with v.net.iso.
It shows distances in meters [400, 800, 1200, and everything further] from a certain point.

My desired output is an Isochrone map. However, I do not know yet how to merge/aggregate isolines into polygons, where each isoline located exactly in a proper polygon with no overlaps etc. Better to say I want to draw an outline from the farthest points of each isoline class. What is the best approach?
I have tried:

extracting nodes from lines Vector > Geometry Tools > Extract nodes, filtering with a certain isoline class, and then Vector > Geometry Tools > Concave Hull.
deploying the QGIS plugin, namely QGIS-ConcaveHull-Plugin.

Both outputs are represented on the image below, numbered appropriately. But there is a problem, nodes from the higher-class isoline group overlap the lower-class isoline group. How can I handle this issue?

I have seen these articles but nothing spectacular:

Drawing outline of point cloud using QGIS?
Creating drive time polygons using Open Source tools?
Computing Walking and Cycling Isochrone maps using QGIS & pgRouting?
How to create contour lines from point data using QGIS?
Create polygon layer from polyline layer


texts: Is it possible to round decimal places for float data in Google Earth Engine?

texts:
I am exporting a large amount of NDVI raster data in Google Earth Engine and was wondering if it was possible to round (my NDVI) float values to 2 decimal places as a way to decrease the size of my exported files (and possibly speed up export time).

I know it is also possible to multiple my data by 100 and convert to integer with round() or int(), then divide by 100 later on if I need the exact NDVI values, but I still wanted to see if I could keep my NDVI values as floats with 2 decimal places to avoid that hassle later on.  Thanks!


texts: Deleting shapefiles and rasters using arcpy

texts:
I am trying to delete a series of files from a gdb and I am encountering suprising difficulties doing so. When I test my procedure in one of the files, let's call it 'yy.shp'This is what I do:

from arcpy import env
from arcpy.sa import *
import os
env.workspace  = env.scratchWorkspace = "path/to/my/GEODB.gdb"
arcpy.Delete_management(r'yy.shp')

And I receive the error:

  Runtime error  Traceback (most recent call last):   File "",
  line 1, in    File "c:\program files
  (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 3995, in
  Delete     raise e ExecuteError: ERROR 999998: Unexpected Error

Notice that if I do:

arcpy.Exists(r'yy.shp')

I get a TRUE result

What am I doing wrong? I have checked and there is definitely no mistake in the path or names of the files. I got these names by dragging and dropping the gdb and file into the python console.


texts: How to query raster in PostGIS within GeoJSON mask only

texts:
I have a raster in PostGIS which I need to query pixel values from, but only within a specific area. So I have a GeoJSON polygon within which I want to extract min/max/avg values of the raster from.

The raster is divided into 640 rid's as well, due to pretiling. How can I combine the values?

So, how can I do this?


texts: QGIS 3.0.1 toggling categories visible/invisible

texts:
If I have a layer with many style categories, all but one or two set to invisible, is there a way to toggle them all visible?

The workaround is to create a duplicate of the layer with everything visible, but is there a better way?


texts:
I working on resource deployment simulation script which works on time enabled features and network analyst.

The resource layer is point feature, each point representing a grouping of units. The count field attribute in this layer should update with time as a unit is deployed (-1) to an incident and returned to home location (+1).

Is it possible dynamically edit field attribute values at specific df time stamps


texts: Exporting NDVI as TIFF file in Google Earth Engine

texts:
I've calculated NDVI using Sentinel data.  Below is the part of my code. This code exporting TIFF file into my google drive, however, when I load this image into ArcMap and trying to convert this raster into shapefile its not happening. This image is opening into ArcMap, but its "Open attribute table" is not enable. I want to share this NDVI layer with other colleagues, so that, they able to calculate area of my NDVI layer. How can I resolve my problem? I'm new to Google Earth Engine
//setting range of my NDVI
var ndviRange1 = 0.344;

// Calculating NDVI

var ndvi1 = img.normalizedDifference(['B8A','B4']).rename('NDVI');
var ndviMask1 = ndvi1.gte(ndviRange1); //geater than

// Updating Mask
var NdviData=ndvi1.updateMask(ndviMask1);

//Exporting this mask
Export.image.toDrive({
      image: NdviData,
      description: 'DGKhanDivCo',
      folder:'DGKhanDivCo',
      scale: 05,
      region: roi,
      fileFormat: 'GeoTIFF',
      maxPixels: 1e12,
    });


texts: Layer layername1 : OGR[3] error 1: REPACK of layer 'layername2' failed. Is there a fix?

texts:
Is there a chance to save my data as they disappeared? Geometry is gone, attribute table is empty. I duplicated layer so probably it is where the problem started.
MacOs Catalina 10.15.7
QGIS3.10
I use QGIS at basic level so is there please anyone who would explain it to me simply?


texts: Looping through folders and select shapefiles by name using pyqgis

texts:
I want to loop through folders and find shapefiles which contains the name "GST" and merging them together using python in QGIS 3.4.3.

In detail:
The folder x contains four folders called VA72, VA73, VA75 and VA76. In each of this folders where a lot other folders (72001, 72002, a.s.o.). In each of them are three other folders, one of them is called "SHP". And in this folder there are several shapefiles. One of the shapefiles is called 72001GST_V2.shp. In every folder (VA72, VA73, VA75 and VA76) the structure is the same.

I want to select all shapefile layers which contain the phrase "**GST_*.shp" and merge them together.

I tried a code I've found after modified it, but after performing it nothing happens.

The code is as follows:

import os
import processing
import glob

Select_folder="C:/*/"
Save_results= "C:/*/merged"
layers = "*GST_*.shp" #select only shapefile layers that contain the phrase GST
os.chdir(Select_folder)

# loop through folders and merge them together
for layers in glob.glob( Select_folder ):
    processing.run("native:mergevectorlayers", {'INPUT':[layers],'CRS':QgsCoordinateReferenceSystem('EPSG:31255'),'OUTPUT':[Save_results+'/'+'GST.shp']})

## load merged GST layer
GST_merged = QgsVectorLayer('C:/*/merged','GST','ogr')
QgsProject.instance().addMapLayer(GST_merged, True)

What is wrong? Has anybody an idea?


texts: Possible to automatically export atlases from QGIS containing PostgreSQL views/tables?

texts:
I've managed to export several atlases automatically from QGIS with a Python and a .bat file, see Export several configured atlases with a Python script

I'm using QGIS 3.10.2 on a Windows 10 PC and I'm connecting to PostgreSQL 9.5.2.

The QGIS project (and layouts/atlases) only contains shapefiles and WMS.

When I try to do the same thing with layouts containing PostgreSQL views and tables, no JPEG files are created. May it be due to PostgreSQL connecting issues?

The username and password for the PostgreSQL/PostGIS connection is stored in the QGIS project (using Authentication manager when creating a new PostGIS connection). (In addition username and password are stored in pgpass.conf.)

I've tried the code from @Matthias Kuhn (@Matthias Kuhn ) (host, tablename, username etc. are not the real ones used in my actual code)

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# linja over bruker vi for å kunne bruke æøå uten advarsler.

import os
from qgis.core import QgsApplication, QgsProject, QgsLayoutExporter, QgsDataSourceUri, QgsCredentials

#from PyQt5.QtSql import QSqlDatabase

# Initialize QGIS Application
QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)

# # Start a QGIS application without GUI
qgs = QgsApplication([], False)
qgs.initQgis()

uri = QgsDataSourceURI()
# assign this information before you query the QgsCredentials data store
uri.setConnection("HOST", "PORT", "DB_NAME", None, None)
connInfo = uri.connectionInfo()

(success, user, passwd) = QgsCredentials.instance().get(connInfo, None, None)

if success:
    uri.setPassword(passwd)
    uri.setUsername(user)
    uri.setDataSource("SCHEMA", "TABLE_NAME", "the_geom")
    QgsVectorLayer(uri.uri(), "<name of view>", "postgres")

    project_path = 'C:\\adhoc\\lst/lstp.qgs'
    output_folder = 'C:\\adhoc\\lst\\n5g\\'
    layout_name = 'n5 NGO48 sone G'

    # Open existing project
    project = QgsProject.instance()
    project.read(project_path)

    print('Project ' + project.fileName() + ' loaded successfully')

    # Open prepared layout that as atlas enabled and set
    layout = project.layoutManager().layoutByName(layout_name)

    # Export atlas
    exporter = QgsLayoutExporter(layout)
    settings = QgsLayoutExporter.ImageExportSettings()
    exporter.exportToImage(layout.atlas(),output_folder, 'jpg', settings)

    # Close the QGIS application
    qgs.exitQgis()

I also had a shot at @J. Montecolo (@J. Montecolo) including connInf, but connection fails.

...
from qgis.core import QgsApplication, QgsProject, QgsLayoutExporter, QgsDataSourceUri, QgsCredentials

from PyQt5.QtSql import QSqlDatabase

# Initialize QGIS Application
QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)

# # Start a QGIS application without GUI
qgs = QgsApplication([], False)
qgs.initQgis()

host = "HOST"
port = 5432
database = "DB_NAME"

db = QSqlDatabase.addDatabase("QPSQL")
db.setHostName(host)
db.setPort(port)
db.setDatabaseName(database)

uri = QgsDataSourceUri()
# assign this information before you query the QgsCredentials data store
uri.setConnection("HOST", "PORT", "DB_NAME", None, None)
connInfo = uri.connectionInfo()

(success, user, passwd) = QgsCredentials.instance().get(connInfo, None, None)
if success:
    db.setUserName(username)
    db.setPassword(password)

if db.open():
    print("Connection successful !")
else:
    print("Connection failed.")
...

Is it possible to automatically export layouts containing PostgreSQL views/tables? If yes, how may I adjust the Python code?


texts: Conversion of raster file into a ASCII grid file using QGIS 3.4.15

texts:
I know this question was already asked on here, but the answer from the link below did not work, at least on this QGIS, as no ASCII option was available when clicking on the files of type option.

QGIS Export Raster to asc


texts: Is it possible to use classes in Google Earth Engine?

texts:
I am specifically asking about programming in Google's Code Editor of the Earth Engine.

I copy-pasted the example class from W3Schools

class Car {
  constructor(brand) {
    this.carname = brand;
  }
}
mycar = new Car("Ford");

and it gives me the following error:

SyntaxError: The keyword 'class' is reserved (2:0)

This is confusing because I would have expected "class" is not defined in this scope., if it would not have supported classes, but the error seems to suggest that class is recognised, but I'm misusing it.

Any suggestions on how to use classes in Google Earth Engine?


texts: Is it possible to link a QGis project with an excel file?

texts:
I have an Excel file with coordinates and I have used this file to plot the coordinates in QGis. I was wondering if it was possible to permanently link the excel file to the QGis project, so that when I add new coordinates to my excel file, these will show up in QGis when I open the project.

Any ideas how to solve my problem?


texts: Calculating convergence angle for point data in QGIS

texts:
I have been a user of ArcGIS Desktop for many years and this function is easily done within ArcMap using the 'Calculate Grid Convergence Angle (Cartography)' tool.

I am wondering if there is a similar function in QGIS.

A little explanation of my process.

I deal with directional surveys for wells on a constant basis. Typically these are reported to the state in their local State Plane coordinate system and as Grid North.

Since i work in a project that contains multiple states however, converting them to True North instead of their typical Grid North (say New Mexico East NAD83 and Texas Central NAD27) is how i keep them accurate in a NAD27 BLM 13N project. The geologics software I use does the proper rotation for them as long as they are true north.

What I have typically done is create a fishnet in ArcGIS and kept the point file that is a byproduct of that.

In that point file i create fields to store the X,Y, and convergence angle for each point, for their local CRS. I then translate each local CRS (keeping the convergence angle of the original CRS) into the combined NAD27 BLM 13N project and can make a grid out of those points for each states wells. I can then sample the wells to a grid and subtract the convergence angle from the azimuth of the wells directional surveys to create True North surveys.

In QGIS I have found how to do the fishnet portion (much easier) however I cant seem to find the same functionality for storing the convergence angle in a field of that shapefile. I have looked in the field calculator under geometry and typing in various searches. I have discovered how to calculate the X and Ys just cant find that darn angle. For NME NAD83 SPCS convergence angles should be around -.7 to +.6 degrees and Texas Central (since Texas went horizontal with its SPCS's) is around -3.5 in the East to +3.25 or so in the West.


texts: Reverse geocode with PostGIS use OpenStreetMap Database

texts:
I want to reverse geocode with PostGIS use OSM database, I use function:

<!-- language: lang-sql -->
CREATE OR REPLACE FUNCTION khoangcach(lat double precision, lon double precision)
RETURNS text AS $BODY$
declare point geometry; rec record; geocode text;
begin point := geomfromtext('POINT('||lat||' '||lon||')', 900913);
select name, ST_Distance(way, point) as dist into rec
from planet_osm_point
order by dist asc;
geocode := rec.name;
return geocode;
end;
$BODY$
LANGUAGE plpgsql;

I think the function is correct but I just have a result for all POINT input.


texts: Leaflet: Resize markers in layer when zoom in

texts:
May I ask how to resize markers in layer when zoom in?

var map = L.map('map', {
     center: [23.364739, 115.187951],
     zoom: 11,
     minZoom: 11,
     maxZoom: 18,
});

var ar_icon = L.Icon.extend({
});

var ar_icon_1 = new ar_icon({iconUrl: '../_image/arrows/line_ar_0.png',     iconSize: [11, 11]});
var ar_icon_2 = new ar_icon({iconUrl: '../_image/arrows/line_ar_1.png',     iconSize: [22, 22]});

var all_testptLayer = new L.GeoJSON.AJAX("../data/test_geojson.php", {
    pointToLayer: function (feature, latlng) {.
        if (feature.properties.num < 0.5)
            return L.marker(latlng, {icon: ar_icon_1});
        else if (feature.properties.num < 1.0)
            return L.marker(latlng, {icon: ar_icon_2});
    }
}

all_testptLayer.addTo(map);

map.on('zoomend', function() {
   var currentZoom = map.getZoom();
   if (currentZoom>12) {
   // what should I do here to access all markers in the all_testptLayer and change the size of all markers to double in size?
   }
});


texts: qgis errpr :The following layers were not correctly generated

texts: I am a novice, is this what is the problem?


texts: SelectLayerByAttribute_management selects all of type of attribute instead of specific one?

texts:
I want to create a tool that allows a user to select all the points within a specific zip code by inputting that zip code. The code below selects ALL of the points with a zip code - which is all of them.

I've tried various iterations of the SQL expression, but this the only one that doesn't cause a crash. Please let me know if you have a solution.

The code below works:

# Author: John K. Tran
# Contact: jtran20@masonlive.gmu.edu
import arcpy
import os
from os import listdir
import csv
print "arcpy imported"

"""
Input CSV can look like:

Lat,Lon,First Name,Last Name
12.34,56.78,Joe,Smith
0.98,7.65,Jane,Doe
65.43,43.21,Bob,Sagat

Make sure 'Lat' and 'Lon' fields are the 1st and 2nd column in the CSV respectively.

"""

incsv = arcpy.GetParameterAsText(0)# Change this to the path of your CSV file.
outfc = r'N:\ReferenceMaterials\ArcGIS Python\Pavement Cores by County\A Interim Shapefiles' + '\\' + arcpy.GetParameterAsText(1) + '.shp'  # Change this to the path of your output FC.

spatialref = arcpy.SpatialReference(4326) # Create the spatial reference object as WGS84. Can modify if desired.

#check for and erase temporary data
testData = r'N:\ReferenceMaterials\ArcGIS Python\Pavement Cores by County\A Interim Shapefiles'
if len(os.listdir(testData)) == 0:
    pass
else:
    for t in listdir(testData):
        if t.endswith('.shp'):
            arcpy.Delete_management(os.path.join(testData, t))

arcpy.CreateFeatureclass_management(os.path.dirname(outfc), os.path.basename(outfc), "POINT", None, None, None, spatialref)

csv.register_dialect("xls", delimiter=",", lineterminator="\n") # Register the dialect for native CSV syntax in Microsoft Excel.
f = open(incsv, "r")
reader = csv.reader(f, dialect = "xls")

headers = reader.next() # Read the first line as the header names.

for header in headers[2:]: # Add fields for remaining columns if needed. Default is TEXT field.
    arcpy.AddField_management(outfc, header, "TEXT")

cursor = arcpy.da.InsertCursor(outfc, ['SHAPE@XY'] + headers[2:]) # Create InsertCursor.
count = 0
for row in reader:
    if count % 1000 == 0:
        print "processing row {0}".format(count)
    Ycoord = row[0] # Make sure 'Lat' is in the 1st column.
    Xcoord = row[1] # Make sure 'Lon' is in the 2nd column.
    newrow = [(float(Xcoord), float(Ycoord))] + row[2:]
    cursor.insertRow(newrow) # Insert point in FC for each row in CSV.
    count += 1

del cursor
f.close()

#set variables for mxd and dataframe
mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, '*')[0]

#import new shapefile into dataframe
newLayer = arcpy.mapping.Layer(outfc)
arcpy.mapping.AddLayer(df, newLayer, 'TOP')

#select attribute (in this case, zip code)
InputFeatureClass = arcpy.GetParameterAsText(1)
InputField = 'zipCode'
InputValue = arcpy.GetParameterAsText(2)
whereclause = """{} = '{}'""".format(arcpy.AddFieldDelimiters(InputFeatureClass, InputField), InputValue)
arcpy.SelectLayerByAttribute_management(InputFeatureClass, "NEW_SELECTION", whereclause)


texts: Building multiple extents as polygons from lat/long list values in R

texts:
I have a loaded .csv table in R containing xlongl, xlongr, ylatb, ylatup values (EPSG:4326). Each row indicates with its combination of lat/long the extent of one rectangled polygon.

> mydata
    xlongl     xlongr      ylatb    ylatup
  1    13.500000  15.500000   6.000000  8.000000
  2     9.000000  22.000000   0.000000 20.000000
  3    18.000000  27.000000  -5.000000  9.000000
  4     9.000000  18.000000  -5.000000  9.000000

Referring to this post entry it is possible to nest the extent function, from the raster library, in as to create a SpatialPolygons object. It seems to be the smarter solution instead of using the more complex way from matrix/data frame(lat/long) -> Polygon -> Polygons -> SpatialPolygon

Could someone help me to develop a coding containing a loop that built the extent as SpatialPolygon object one row after another? And that can finally merge the created single SpatialPolygon objects to one by maintaining the original order of the table?


texts: How to add offset to layer in OpenLayers?

texts:
As the title suggests,I have GeoJSON vectors on OpenLayers map where user should add X/Y offset (i.e 5 meters north,5 meters east).I plan to add some control for it.

Since vectors are drawn map with OpenLayers.Rule object,as I searched for some offset option,then all I found was graphicXOffset & graphicYOffset which eludes normal vectors but only used for external graphics.

Is there any easy way to accomplish this feature on client ?


texts:
Whenever I attempt to run use run the Voronoi Map in the explore data option of the Geostatistical Analyst tool an error message as shown below pops up.

When I run the debugger all the errors say is "Cannot find or open the PDB file".
I am not sure how to address this issue.  I have already tried renaming the normal template to OLDnormal for ArcCatalog and ArcMap and have reset my TEMP folder.  Has anyone else had any success in addressing this error?


texts: I'm trying to join an excel file to ArcMap. The information in the excel file is clean and doesn't contain any characters that might hinder the joining process. However, when I join the file, all of the information in the attribute table comes up as "null", and the information in the top row of the attribute table is the first county in the first cell of my excel file. I was wondering if there was a way to fix this.


texts: Registered view gets buggy when database link added

texts:
I have a view in an Oracle 12c geodatabse:

create or replace view parcel_vw as (
select
    p.*
from
    parcel p
left join
    for_gis.gis_locations@gistomax loc  --<-- Join to database link
    on p.propertyid = loc.location
)

[there are 680 records in the view/fc]

The view has a database link to a table in a separate Oracle database.

Note: The join to the dblink does not propagate-out more records; the relationship is 1 to 1.

This is what happens when I use the view in Arcmap 10.7.1:

Registered=n; dblink=y; WORKS
Registered=y; dblink=n; WORKS
Registered=y; dblink=y; BUGGY

In #3, when I try to select all features, I get buggy results:

Occasionally: ArcMap crashes
Consistently: Not all of the records can be selected with the select tool in the map.
Consistently: None of the records can be selected in the attribute table with the Select All tool.

Why does the view (with dblink) become buggy when registered with the geodatabase?


texts:
I have a problem with the Profile Tool plug In.
I have a layer which contains the X,Y and elevation of a cross section.
I add new points at both extremities, entering only the elevation value.
When I use the Profile Tool, it plots me the profile between my first points and then the profile between both extremities.  Which gives a 'mirror' effect (see below). As if there was an order for plotting the profiles.
Is it possible to plot the profile of all the point from one extremities to another ?


texts: Calculate land use percentage within an area?

texts:
I'm interested in calculating land use percentage for three layers (forest, agriculture and urban) within the context of a water basin area.

I've managed to merge the land use layers with Union tool and then clipped to the area of the basin. But having trouble thinking how to do it.

The result I'm thinking I want is an image showing the area with the three different land use polygon visible and then percentage in the attribute table.


texts: What layers are produced by SAGA Potential Incoming Solar Radiation?

texts:
I estimated the Potential Incoming Solar Radiation in SAGA. The tool creates a file "direct insolation" and "diffuse insolation". When I export these files as .tif, a layerstack of 3 raster layers is produced for each insolation file.

In the documentation of the tool I can not find an advice, what the layers in particular mean. Does anyone knows how to interpret those layers?

Here is a screenshot of the layers from direct insolation


texts: Is it acceptable to assume there will be no 'walk' in coordinates transformed back and forth between two datums?

texts: Given that you use the same datum transformation, you should end up with the same coordinate regardless of how many times you convert between to datums.  Is this a correct assumption?  Is this applicable for both equation-based and grid-based transformations?  For example: Convert NAD27 to NAD83 via NADCON then from those NAD83 to NAD27 via NADCON, and again, etc.


texts: Implementing non-default keyboard shortcuts in Openlayers2

texts:
Following on from this thread, I am trying to implement @Aragon's suggestion and add some keyboard shortcuts to OpenLayers via a handleKeyEvent method in OpenLayers.Handler.Keyboard.  This ought to be a simple task that many people implement but for some reason I'm just not getting it.

var keyboardControl = new OpenLayers.Control();

var keyHandler = new OpenLayers.Handler.Keyboard(keyboardControl, {
    handleKeyEvent: function(evt) {
            alert("You pressed a key!");
            switch(evt.keyCode) {
                case 87:
                    doSomething([0, 1]);
                case 83:
                    doSomething([0,-1]);
                case 65:
                    doSomething([1,0]);
                case 68:
                    doSomething([-1,0])
            }
        },
});

function doSomething(vals){
    alert("Something will happen here eventually...");
};

map.addControl(keyboardControl);

I get no errors but neither of my alerts get fired.  What am I doing wrong?


texts:
I am a beginner to this field and I am using Open Data Cube for the first time and after following the installation steps from (https://datacube-core.readthedocs.io/en/latest/ops/db_setup.html)
I am getting the error No datacube config found when I run datacube -v system init.


texts: Getting bounding box from Point type GeoDataFrame to getting data from Osmnx

texts:
I've a point geodataframe which has only point features inside of it.

fishnet_Id           ...                          geometry
0           0           ...            POINT (403486 4540640)
1           0           ...            POINT (403736 4540640)

I want to create a bounding box polygon feature from the easternmost, southernmost, northernmost and westernmost point and then pass it into osmnx.gdf_from_bbox method.

I could do it via finding the biggest/smallest lat and biggest/smallest lon values and then make a column from these values and pass these values into the osmnx function manually but I want to do it more code-efficient way.

(Right now its look like this)

import geopandas as gpd
import pandas as pd
import osmnx as ox
from fiona.crs import from_epsg

origins = gpd.read_file(r"data\Origin_Clip.shp")
destinations = gpd.read_file(r"data\Dest_Clip.shp")

origins["geometry"] = origins["geometry"].to_crs(epsg=4326)
destinations["geometry"] = destinations["geometry"].to_crs(epsg=4326)

origins["lat"] = origins["geometry"].x
origins["lon"] = origins["geometry"].y
destinations["lat"] = destinations["geometry"].x
destinations["lon"] = destinations["geometry"].y

points = gpd.GeoDataFrame(pd.concat([origins, destinations], ignore_index=True))

xmin, xmax = points["lat"].min(), points["lat"].max()
ymin, ymax = points["lon"].min(), points["lon"].max()

G = ox.graph_from_bbox(ymax, ymin, xmax, xmin, clean_periphery=False)


texts: Generate dynamic input with python script

texts:
I'm using 10.2.2 ArcGIS, I want to create a python script and have the following script:

I'm new to this topic
as I can use sys.argv [] or GetParameterAsText?
for can I make dynamic input coordinates and the name of the output shapefile

arcpy.env.workspace = 'C:\Users\Modelo IDEAL\Puntos.gdb'
arcpy.env.overwriteOutput = True
sr = arcpy.SpatialReference('WGS 1984')
arcpy.env.outputCoordinateSystem = sr
latitude = 14.024982
longitude = -89.962368
p = arcpy.Point(longitude, latitude)
pg = arcpy.PointGeometry(p)
arcpy.CopyFeatures_management(pg, 'salida.shp')


texts: get rgb images from sentinel-2 using gdal

texts:
i have some troubles to get RGB true colors image from sentinel satellite bands

wich produce this

to me it looks not natural at all

also i have tried gdal to do the merge with :

gdal_merge -o merged.tif -of GTIFF B04.jp2 B03.jp2 B02.jp2

but gave me an unusable file, at least in microDEM .

so i am wondering what's wrong here ?


texts: Building a Polygon from an array of points in Arcobjects and C#

texts:
I have seen this code in Arcpy which is used to create a Polygon from a string of Points(X , Y).Can you please let me know how I can us something like this which can generate polygon from an array of X,y values in C# and Arcobjects?

I already take alook at Arcobjects 10 .net framework sample code sample of Building a Polygon Using Points . but it was too confusing for me.

I really appreciate your time if you can point me to a tutorial or a sample which contains a process of polygon creation same as what was in arcpy code.


texts: Does sum line length measures geodetic or planar distance?

texts:
In QGIS i created a grid (via "create grid"), which i later split into layers of horizontal and vertical lines. Afterwards i calculated the respective cumulative line length for polygons using "Sume line length". The results i achieved differ especially for horizontal lines quite strongly (~ 25% less) from expected values (in the case where i could compare, based on another publication (values for China)).

One suggestion was, that the difference may be due to measuring planar vs. geodetic distance. Does anyone know how QGIS measures distance with this tool? May the result be subject to the projection I use?


texts:
I am trying to import a geojson file with the US States.
The geojson file has 3 fields. gid,state_abbr, and geometry.

I submit this code

bin/geomesa-cassandra ingest --contact-point 127.0.0.1:9042 --key-space mykeyspace --catalog mycatalog --converter state-json --spec state-json states.geojson

and get back

INFO  Creating schema state-json
INFO  Running ingestion in local mode
INFO  Ingesting 1 file with 1 thread
[                                                            ]   0% complete
0 ingested 57 failed in 00:00:01
INFO  Local ingestion complete in 00:00:01
INFO  Ingested 0 features and failed to ingest 57 features.

My converter file is

"state-json" = {
  type         = "json"
  id-field     = "$gid"
  feature-path = "$.features"
  fields = [
    { name = "gid",   json-type = "integer",  path = "$.gid",             }
    { name = "state_abbr", json-type = "string",   path = "$.state_abbr", }
    { name = "geom",  json-type = "geometry", path = "$.geometry",    }
  ]
}

and my spec file is

"state-json" = {
  attributes = [
    { name = "gid", type = "Integer"    }
    { name = "state_abbr",  type = "String" }
    { name = "geom",   type = "Geometry"  }
  ]
}

Is there anything I am doing wrong?


texts: Symbolizing Certain Attributes in a Polygon and the Excluding One in QGIS

texts:
I have a map full of 78,000+ parcels. I want to fill in the color of all parcels owned by the US Government. One problem is they are named several different ways in my data ("United States of America", "US Forest Service", "United States Forest Service", etc.). So I built a formula that looks like this:
"ACCT_NAME" = 'UNITED STATES FOREST SERVICE' OR  "ACCT_NAME" = 'U S NATIONAL FOREST' OR "ACCT_NAME" = 'UNITED STATES OF AMERICA'

This works and narrows down all the entries to 66.
Now out of 66 government owned properties, I have one of the parcels that is owned by "UNITED STATES FOREST SERVICE" that I don't want to display. So I tried to eliminate this single parcel by another attribute, so I added the following:
"ACCT_NAME" = 'UNITED STATES FOREST SERVICE' OR  "ACCT_NAME" = 'U S NATIONAL FOREST' OR "ACCT_NAME" = 'UNITED STATES OF AMERICA' OR "PIN" != '7626736067'

This doesn't work. I tried adding "IS NOT" "<>" and "!=" as the operators, but none seem to work.
Any suggestions?


texts: Making polygon layers always have area fields in QGIS

texts:
I am making Delaunay triangles of several (over 40) combinations of points and it is very tedious to output new shapefiles for each combination without even having the area of the triangles by default, which is what I am most interested in.
Is there a way to make an area field by default for all new polygons and besides can I use tools like Delaunay triangulation to add new polygons to one shapefile without having to make new shapefiles each time?


texts: How to fill unconnected areas in polygon with the same attribute

texts: I am trying to fill in the gaps with the same attribute.  These arrowed areas belong to a city called 'Westbank', but the thing is, these areas are unconnected, how can I make a polygon out of areas by making sure these 2 areas will be unconnected? Can I even do this in QGIS3?


texts: Layer property refuses to act as anything other than integer when named "SERVICE"

texts:
I've been trying to add a new property to a vector layer named SERVICE, this property should be a string. I am able to create this property no problem, and the type is stated to be QString in layer properties.

However, when I then open the feature attributes window to edit this property on a feature, the field appears and acts as an Integer, accepting only numbers and showing the adjustment arrows.

This strange behaviour disappears if I change the name of the attribute, even if I change it simply to "Service". For some reason the all caps version of this name seems to cause the field to act as an integer, even though it remains registered as a string in the layer properties.
I'd be interested to see if others are able to reproduce this behaviour, and if anyone knows why it is the case.


texts: How to "save as" a layer creating new data files, with QGIS 3.4?

texts:
With QGIS 2 you could right-click on a layer; choose "save as"; and create a new layer based on new data files that would thus be created (as contrasted with "Duplicate", which would create a new layer based on the same old data files.)

With QGIS 3.4.4 I see "Copy Layer", which I can then "Paste", but that seems to create a new layer with the same name and based on the same old data files.

QGIS 3 also offers "Duplicate Layer", but that also creates a new layer based on the same old data files.

So how do we now create a new layer based on a new set of data files?


texts: Automatic selection of colored regions when digitzing choropleth maps

texts: I am new to working with GIS. I have a raster choropleth map similar to the one displayed below. I want to create a shapefile from the map, where each layer contains all the countries that have the same color. In this example, Brazil and Argentina would be in the same layer. Is there a way to automatically select all the regions that have the same color? Alternatively, can I select an area that is enclosed by the black lines, e.g. Brazil? In photoshop, there is a "quick-selection" tool that does this. I am using QGIS.


texts: Polygon which is drawn only with 2 points is not saving in ArcGIS Online service

texts:
I have a polygon layer in ArcGIS Online,

When I draw a polygon by taking P1 as starting point and double click on the P2 so polygon draw end event will come from draw tool.

The below is the JSON format which I used to addFreature.

[
{
  "attributes" : {
    "attr1":"10008201",
    "attr2":"9",
    "attr3":"PMS"
  },
  "geometry" : {
    "type":"polygon",
    "rings":[[[8487016.789333614,1464112.0521501263],[8668631.168539183,1622489.5747570037],[8668631.168539183,1622489.5747570037],[8487016.789333614,1464112.0521501263]]],
    "_ring":0,
    "spatialReference":{"wkid":102100,"latestWkid":3857}

  }
}

]

The addFeature service says sucess. But the geometry of that feature is not saved in the ArcGIS Online layer.

So, I am thinking the polygon geometry is invalid so its not saving.
But how to validate it from JavaScript API


texts: GRASS 6.4.3 Crashes on Startup; libintl-8.dll cited

texts:
I recently installed QGIS 2. I used the advanced OSGeo4W installer that I downloaded from the QGIS download page. GRASS is included in the stack. When I attempt to open GRASS, I receive the following error message:
The procedure entry point libintl_setlocale could not be located in the dynamic link library libintl-8.dll.
I checked the WinGRASS errors page on the GRASS wiki and did not discover any precedent for this problem either there or by a Google search.
My operating system is Windows Vista. I've installed all recommended packages as outlined on the GRASS wiki.
Has anyone else encountered this problem? And how might I go about resolving it?
Here's a running list of potential solutions that have not corrected the problem. I'll update as I get the time to try them:

Run Microsoft .NET cleanup utility and reinstall .NET 3.5.
Run system file checker from command prompt (>sfc /scannow).


texts:
I generated, from point data, a kernel density map using GRASS, and I would like to identify the 95% volume contour in it.

Is it possible in GRASS?


texts: Execute JavaScript code in popup for Mapbox?

texts:
I'm trying to embed tweets into the popups for a map I've made in Mapbox. The text of the tweets is displaying properly, but not the styling — the embed code contains a script tag to load a widget for styling, and I would guess that this is not being executed by the "sethtml" function. How can I execute the JavaScript when the popup is opened in Mapbox?

I've run into this problem in Leaflet before, and solved it with the following function.

map.on('popupopen', function(e) {
    $.getScript("https://platform.twitter.com/widgets.js");
});

I admit I'm new to Mapbox, but I thought this script would also work for mapbox GL JS. I have the following (not working) code for GL JS:

// add markers to map
geojson.features.forEach(function(marker) {

  // create a HTML element for each feature
  var el = document.createElement('div');
  el.className = 'marker';

//get JS widget and move marker down
  map.on('popupopen', function(e) {
    $.getScript("https://platform.twitter.com/widgets.js");
    var px = map.project(e.popup._latlng);
    px.y -= e.popup._container.clientHeight/1.5;
    map.panTo(map.unproject(px),{animate: true});
});

  // make a marker for each feature and add to the map
  new mapboxgl.Marker(el)
    .setLngLat(marker.geometry.coordinates)
    .setPopup(new mapboxgl.Popup({ offset: 5 }) // add popups
    .setHTML(marker.properties.Embed + marker.properties.Name))
    .addTo(map);

});


texts:
I have about 50 pdf maps that I would like to batch georeference and/or georectify. I've saved one as a jpg and georeferenced it manually and would like to automate the rest which are in the same location and have the same extent.

The density gradient will change in each case, but the map (saved as jpg now) will have the same extent for all 50 cases.

Does anyone know of a link to a step-by-step tutorial that deals with batch processing in Arcmap 10.2 only? I'm thinking of using the Warp from File code or possibly the Autoregistration option, but am not sure where to start?

The sample code is:

import arcpy
arcpy.WarpFromFile_management(
"\\cpu\data\raster.img", "\\cpu\data\warp_out.tif",
"\\cpu\data\gcpfile.txt", "POLYORDER2", "BILINEAR")


texts: Update one table from another based on spatial join

texts:
In beginning i had two tables of points and i wanted create from one of this POINT table new table CIRCLES which will contain circle polygons 300 metres around points. This step i have done here:

 INSERT INTO circles (geom) (
   SELECT ST_Buffer(point, 300, 'quad_segs=8')
   FROM points
);

Now I have two tables. First table CIRCLES contain 10 polygons (circles) and each circle have unique serial number (1,2,3,..10). Second table POINTS contain 500 points. How can I, get this result: If any point from table POINTS is located inside polygon (circle) from table CIRCLES, so in table POINTS in column PLGN we will write value of serial number of circle in which is this point inside.

In this task we don't consider case, that one point can be inside two or more circles.


texts:
I can't remove item from the legend (lyr called "atikot") using arcpy with this code:

import arcpy
from arcpy import env

removeLyr = arcpy.mapping.ListLayers(mxd, "atikot")[0]
mxd = arcpy.mapping.MapDocument(r"C:\Project\Project.mxd")
legend = arcpy.mapping.ListLayoutElements(mxd,"LEGEND_ELEMENT")[0]
for lyr in legend.listLegendItemLayers():
    if lyr.name == "atikot":
        legend.removeItem(lyr)
del mxd

The code does not work. Does anyone have an idea as to what's wrong with the code?


texts:
I wish to apply my action on pyqgis (in a plugin).

In my QGIS plugin, I created an action that I apply to a certain layer. The action appears well on my layer but it is not automatically applied, I have to do it manually.

I think I should use the run() function of the QgsAction class.

layer_hypothese = iface.activeLayer()
action_layer_hypothese = layer_hypothese.actions()
actionDescription = 'Graph radar '
actionCode = 'from qgis.utils import iface; iface.messageBar().pushMessage("Feature ID: [% $id %]", 0)'
scope = {'Feature','Layer','Canvas','Field'}

action_graphique = QgsAction(QgsAction.GenericPython,actionDescription,actionCode,'',False,'',scope)
action_layer_hypothese.addAction(action_graphique)


texts: Selected features with a join are identified improperly from a service

texts:
We couldn’t figure out why a service fails to provide correct identification for the selected feature. However, the scenarios below may highlight the issue:

In the “UrbanMasterPlans_WithJoin.mxd” file, the attribute table of the “UrbanMasterPlans_PNA_Approved” layer is joined with the “Regulations” table, then the ArcGIS Server fails to identify all the selected features.

In the “UrbanMasterPlans_WithoutJoin.mxd” file, the attribute table of the “UrbanMasterPlans_PNA_Approved” layer is not joined with the “Regulations” table, then the ArcGIS Server is able to identify all the selected features.

Here is data link to test the problem:
http://www.mediafire.com/download/ev5wthdiijhp51n/PublishingIssue.zip

What might be the issue here? Why does the join affect the Identify results when applied in a service?


texts: ST_Contains when linestring crossing itself

texts:
I have a linestring which crosses itself, like in the picture below. There isn't any vertice at the crossing.

It is a MultiLinestring, so when I want to get the startpoint, the endpoint, or the azimuth (for instance), I get an empty result. I tried to use ST_Linemerge, but for the same result. I think the problem come from the fact that this linestring is crossing itself since I don't have that problem with the linestrings (in the same table) which don't cross themselves.

I am trying to get xmin, xmax, azimuth... in order to establish a process for an atlas in QGIS. I would prefer to manage that in PostGIS, do you have any idea on the way to go?

In my experience these self-crossing linestrings are a big source of problems (with ST_Dump for example), problems that nor me nor my colleagues ever managed to solve, so any insight would be very welcome.