texts:
I'm new to GIS and having a basic beginner problem.  I've got data in the form of a geopackage .gpkg from GADM.org.  I can read in and use the data without problem using Geopandas (via built-in Fiona).

I've added and deleted some columns to the data, and changed some column names to ones more intuitive to me (for ease of later use and consistency across files).  Minor edits.

Now I want to export the data to a file format that is easy to later read into Python as a GeoPandas dataframe.  However, whenever I try

CityShapeFile.to_file('CityShapeFile_edited.gpkg',encoding='utf-8')

No matter which format I choose (e.g. gpkg, json, shp, ...) I don't get a single file, but rather a folder with 5 files: *.cpg, *.dbf, *.prj, *.shp, and *.shx.  I'm sure those files together contain all the data I need, but in order to read the data back into Python I really need it in one file (as far as I know).  Considering the original data I downloaded is in a single gpkg file, this seems possible and probably normal, but how to do it?

I read in the Fiona manual that it can write zipped shapefiles, but I couldn't find any simple example of doing that with a GeoPandas dataframe, nor am I sure whether that can be read in correctly. So,

  How do I write a GeoPandas dataframe into a single file (preferably JSON or GeoPackage)?


texts: Create columns for specific "other tags" of OSM PBF files in QGIS

texts:
When downloading PBF files of Openstreetmap via Geofabrik, the file contains a field called "other_tags" which includes all the OSM tags for this feature. It is a String field and looks like "cycleway"=>"track","lit"=>"yes","maxspeed"=>"30","sidewalk"=>"both","source:maxspeed"=>"DE:zone30","surface"=>"asphalt".

I now want to create new columns for only some specific tags. Lets say "cycleway" and "maxspeed". I need to do this in QGIS.


texts: Restricting information to metadata on GeoNetwork when linking online resource

texts:
How to hide information to the public users when linking online resource using a file on a local network?  The official Geonetwork documentation has a section on the topic here, but the section is empty even though the section Linking website, web services, … using URL referes to it.
The question was asked a year ago but was never answered.


texts: Layer is being add to 2 groups instead of 1

texts:
I'm trying to create 3 group layers, depending on parameters of two lists, but reg_imobiliarios layer is added to auxiliar and referencia groups.

In short, all the layers are in the correct group, except reg_imobiliarios, that I would like to be only in referencia group.

Note: I'm using postgis tag, because maybe is a problem in my database, but I think is improbable.

This is the legend before:

And after running the code:

The results of Python console:

editableLayers = {'aux_area':'auxiliar', 'aux_invasao':'auxiliar',
 'aux_linha':'auxiliar', 'aux_ponto':'auxiliar', 'aux_unificar':'auxiliar',
 'contratos':'ctgeo', 'cotas':'ctgeo', 'imoveis_alienados':'ctgeo',
 'laaa':'ctgeo','lote_legal':'ctgeo', 'cotas_reg_imobiliarios':'referencia',
 'reg_imobiliarios':'referencia'}
deleted_layers = {'aux_area': 'auxiliar', 'cotas': 'ctgeo', 'reg_imobiliarios': 'referencia'}
root = QgsProject.instance().layerTreeRoot()
auxiliar = root.addGroup('auxiliar')
ctgeo = root.addGroup('ctgeo')
referencia = root.addGroup('referencia')

for layerTable, layerSchema in editableLayers.items():
    for tableDeleted, schemaDeleted in deleted_layers.items():
        if layerTable == tableDeleted:
            print(layerTable)
            print(tableDeleted)
            while True:
                print('entrou no while')
                layer = QgsProject.instance().mapLayersByName(tableDeleted)[0]
                QgsProject.instance().removeMapLayer(layer.id())
                uri.setConnection("LPOSTDES-CL2", "5432", "ugpi", user, passwd)
                uri.setDataSource(schemaDeleted, tableDeleted, "geom")
                layer = QgsVectorLayer(uri.uri(False), tableDeleted, "postgres")
                QgsProject.instance().addMapLayer(layer)

                if layerTable == layer.name() and layerSchema == auxiliar.name():
                    print('entrou no auxiliar')
                    print(layerTable)
                    print(layerSchema)
                    vlayer = root.findLayer(layer.id())
                    clone = vlayer.clone()
                    auxiliar.insertChildNode(0, clone)
                    root.removeChildNode(vlayer)
                elif layerTable == layer.name() and layerSchema == referencia.name():
                    print('entrou no referencia')
                    print(layerTable)
                    print(layerSchema)
                    vlayer = root.findLayer(layer.id())
                    clone = vlayer.clone()
                    referencia.insertChildNode(0, clone)
                    root.removeChildNode(vlayer)
                elif layerTable == layer.name() and layerSchema == ctgeo.name():
                    print('entrou no ctgeo')
                    print(layerTable)
                    print(layerSchema)
                    vlayer = root.findLayer(layer.id())
                    clone = vlayer.clone()
                    ctgeo.insertChildNode(0, clone)
                    root.removeChildNode(vlayer)

# Test to confirm if the layer is editable (connection is ok)
                layer.startEditing()
                if layer.isEditable():
                    print('entrou no teste do lapis')
                    iface.vectorLayerTools().stopEditing(layer,False)
                    break


texts: I'm working with a tiff image in arcgis projected in UTM 33N and I need to obtain pixel coordinates of a specific point in the image with known metric coords (UTM 33N). How can I get them?


texts:
To run a simple interest query in Mapinfo I would use

Select * from _1_10000_map_sheet_key_map, Carpentaria_Shire where _1_10000_map_sheet_key_map.Obj Intersects Carpentaria_Shire.Obj into Selection

In QGIS in Add--Virtual Layer and in DB Manager I keep getting errors

select * from _1_10000_map_sheet_key_map where st_intersects (_1_10000_map_sheet_key_map.geometry, Carpentaria_Shire.geometry) as Selection

gives the following error

  Query execution error on CREATE TEMP VIEW _tview AS select * from
  _1_10000_map_sheet_key_map where st_intersects (_1_10000_map_sheet_key_map.geometry, Carpentaria_Shire.geometry) as
  Selection: 1 - near "as": syntax error

The same from Virtual Layers and trying to use geom_from_wkt I get

  Query preparation error on PRAGMA table_info(_tview): no such column:
  Carpentaria_Shire.geom_from_wkt

What's the expected structure of the query in QGIS SQL implementation?

Additionally can this be done in the Advanced Filter expression?


texts:
I'm working with Quantum GIS 1.8 Desktop on Ubuntu 11.10. Two days ago the program crashed while I was working with some WMS layers and then it hasn't been working any more. When I open Quantum GIS Desktop the program starts, the logo appears but nothing more happens and then Ubuntu suddenly quit the application.

I have tried several times both to remove Quantum GIS and to reinstall it, but the problem remains.

Could you help me?


texts: Creating Polys and Displaying geo points within. ( Where to start? )

texts:
Hey all,
  I'm somewhat of a newb when it comes to creating mapping solutions. I've stayed on a pretty high level playing with openlayers and pulling data from various webservices.

The first part of my question is the realm of how do I compute polys?

I understand poly basics, such that it's a 2 dimensional shape that you can create a boundary with. For example, when you enter your zip code in the following example. The red line illustrates the boundaries, which is understood as a poly, correct? Could you link me information on this concept, and how to create polygons and how to work with them. I'm very fluent in python, and I think there is a GIS-python module/postgre-gis?

The second part is, how do I go about matching coordinates in said polys? What is the algorithm I need to make this happen? Is there an open-source project that I can use to calculate this? Preferably in c++ or python.

Third, where can I obtain datasets of Country, State/Territory, County/Providence, City/Town, Neighborhood.

Thanks in advance!


texts: Changing coordinate system says that geotransformation is invalid?

texts:
I am planning on using the MODIS burned area products which were produced by Maryland University. These come as either polygons or rasters. In the user guide  it says

  A shapefile version of the MCD45 product is derived from the monthly Geotiff version by the  University  of  Maryland.  The shapefiles  are  available  with  the  same  projection (Plate-Carrée) and geographic extent as the Geotiff sub-continental windows

The rest of my data are rasters (WGS 1984 UTM ZONE 38N).

When I load the MODIS fire monthly rasters into ArcGIS it says that the system coordinate system is different.

However, when I tried change the coordinate system the results output say that the geotransformation are invalid.

What am I doing wrong?


texts: gvsig raster image misplaced in leaflet, projection problem?

texts:
I am trying to show in leaflet a geotiff raster image that was generated in gvSig.
The basemap I am using in leaflet is a wms server with projection EPSG:3857.

Leaflet needs to know the upper right and bottom left corners in latitude and longitude coordinates which I get from the tfw file.

When I add my raster image the south west and north east corners are displayed in the correct position (the ones I specify to leaflet), but the other two are not. And most of the image is not displayed at the correct position either.
It seems there is a problem with the projections.

gdalinfo of the raster file shows:

Driver: GTiff/GeoTIFF
Files: ModCombustible.tif
Size is 600, 945
Coordinate System is:
PROJCS["ETRS89 / UTM zone 30N",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","25830"]]
Origin = (713179.879402380902320,4381253.375599999912083)
Pixel Size = (5.001750158730428,-5.001750158730441)

I tried to use gdalwarp to reproject the raster in EPSG:3857 with the command:
 gdalwarp -t_srs EPSG:3857 ModCombustible.tif ModCombustible_3857.tif

but in the tfw file I am getting negative UTM easting coordinates which doesn't make any sense to me.

What should I do to properly show the raster image in leaflet?

EDIT: The negative easting coordinates are correct for the EPSG:3857 if they are west of London. Using the link provided by the user30184 I could translate the coordinates to latitude and longitude and show the image at the correct position.

EDIT2: To clarify basic concepts which were not clear to me: as IvanSanchez mentioned in the comments, UTM (Universal Transverse Mercator) is a system of 60 different projections using false easting and northing. In the other hand Web Mercator (EPSG:3857) is not UTM but uses easting and northing coordinates to define a position. It is the de facto standard for web map applications.


texts: How to display image from SQLite database at each point? (QGIS)

texts:
How do I display images at specific points, through the raster image marker symbology, from a standalone BLOB attribute table (SQLite) loaded into QGIS?
To briefly summarise: I have loaded a SQLite database and a geopackage database into QGIS. This geopackage database contains points spread over all of Denmark, containing information like GNSS station names, gps nr, coordinates, height measurements etc.
Then, as a supplement to these points, I have this SQLite database that contains images tied to each point. The SQLite database loads in fine, and can join layers with the geopackage database. But I can't seem to figure out how to display the images tied to each point.
So far, I've managed to display a custom image that I loaded down from Google at each point, through the raster image marker, roughly following this guide:
Guide to adding pictures through raster image marker
But I want to load the images from the database (geo_ref skitse) that has unique images tied to each point instead, which is a bit more tricky.
I have made 2 screendumps of how my situation looks as of now, if that is of any help.


texts:
I am drawing points that have different tags, for example, one is a Manhole other is a Headwall. I have my symbology setup in such a way that point displays differently based on tags. Now as I am digitizing hundreds of features I dont want to tag each one of them after drawing, I want to draw the features with a default tag.

In ArcMap this can be done with feature templates. See below:

Is there a similar function in QGIS?


texts: Getting exact maximum daily temperature using remote sensing?

texts:
I have a problem I want to find the land surface temperature of Bangkok using landsat8.

I have everything to get it but the problem is landsat images are obtained during morning time but the maximum temperature will be afternoon period.

How can I get the exact maximum temperatures of a day?


texts:
I am trying to force a tab from a tab container to open to a specific tab once a new MXD is opened. Currently it will stay on a previously selected tab when I use the catalog(docked) to open a new MXD. My question is how do I tell if a new MXD is loaded?

Thanks
-Chris


texts: QGIS Spatialite import fails when data is projected

texts:
I'm trying to import a shapefile into a spatialite database via drag & drop. This works fine while the data is in a geographic coordinate system, but I need to do some measurements and therefore need projected data. However, I get an error when I try to import the projected data. The error message isn't very helpful either. How can I solve this?

I also tried the DB manager and geht this error:

Spatialite-GUI: no "Map Preview" option

EPSG 54004:

EPSG 3857:


texts: Calculating Image boundary / footprint of satellite images using open source tools?

texts:
I need to create the polygon outlines of several single layer raster images, not the extent/bounding box, but the area without the nodata values, as shown here: Creating shapefile showing footprints of Rasters?.

In the answer to the question above, the Image Boundary plugin was mentioned, however, I don't find it in the QGIS 1.8.0 Lisboa, running on Ubuntu.

Is the tool still available?

If not, is there a way to do this using opensource tools, either gdal, R, QGIS, GRASS, or similar, not with ArcMap?


texts: Calculating median of multiple polylines in ArcGIS for Desktop?

texts:
I am working with several polylines of mean monthly Arctic sea ice extent for the month of September; one line per year from 1979 - 2010.

I would like to summarize the polylines by showing the median (of the mean monthly) ice extents over two different time periods (1979-1999 and 2000-2010).

Is there at tool that will calculate this new median line for me?


texts: Unable to view data in Mosaic Dataset in full extent?

texts:
I have created a Mosaic Dataset containing Raster Datasets which point to TIF files shared at a remote Windows server.

I can see that the Dataset created is fine as I could see grid cells but I could not see the actual Image from Full extent. I have to zoom it to see the images.


texts: How to save changes in Feature-Geometry with GDAL/OGR and Python? (Segmentation fault)

texts:
I have a problem with my code:

driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open('merge.shp', 0)
if dataSource is None:
    print 'Could not open file'
    sys.exit(1)
layer = dataSource.GetLayer() #layer = Shapefile

numFeatures = layer.GetFeatureCount()
i = 0
intersec_count=0
while i < (numFeatures-2):
    feature1 = layer.GetFeature(i)
    geom1 = feature1.GetGeometryRef()
    j = i+1
    while j < numFeatures:
        feature2 = layer.GetFeature(j)
        geom2 = feature2.GetGeometryRef()
           if geom1.Distance(geom2)<0.0001 and geom1.Distance(geom2)!=0:
              buff = geom1.Buffer(0.0001)
              linkline = geom2.Intersection(buff)
              linkpoint = linkline.GetPoint(0)
              geom1.AddPoint(linkpoint[0],linkpoint[1])
              feature1.SetGeometry(geom1)
        j = j+1
    i = i+1
dataSource.SyncToDisk()
dataSource = None

When running this script the folloging Message appears:
"Segmentation fault (core dumped)"

All I'm trying to do, is to add a Point (linkpoint) to a LineString (geom1) and save these changes.

I hope someone can help me :)


texts: I am using QGIS 3.4 Madeira. I want to use difference tool from SAGA from the processing tool box since my data has some invalid geometries and want to fix the issue later after using difference tool so that some invalid geometries are reduced. However, I can't find the tool while searching. I could find the tool while using QGIS 2.18 Las Palmas. Does any body have idea how can I find find SAGA difference tool in QGIS 3.4?


texts: PostgreSQL/PostGIS: Adding attributes to a table that I can use to query with

texts: Currently I have a PostGIS database that is taking in shapefiles. However, these shapefiles are converted into tables when they are imported into the database. Each shapefile deals with a farm name, field name, and year. Is there a way that I can add these attributes to be associated with the table so I could query my table by farm name and field name?


texts: Escaping cut of western/eastern hemispheres in QGIS?

texts:
I am very new to QGIS and GIS in general, so I am struggling with the following: when I do a map of Russia, I get the part next to Alaska being drawn far left from Europe, as it lies behind 180o... I could not find any CRS that depicts Russia "in full"...

Can someone please explain to me how wrap the map in being "one chunk", not being split by the 180 longtitude?


texts: Can I use GeoJSON GeometryCollection in QGIS?

texts:
I want to create a GeoJSON that contains multiple geometries to be used in QGis later on.
Now I red through the specs here: http://geojson.org/geojson-spec.html#geometrycollection

When I simply copy and paste the GeometryCollection example into a file and add it to QGis nothing shows up.

Am I doing something wrong or is it a QGis Bug?


texts: I have a GEOTiff file that has elevation data associated with it .  It seems like the resulting png/jpg files I get after feeding it to gdal2tiles does not have elevation data. Same thing with the tms tiles provided by Geoserver.  Does tiling a GEOTiff file strip away the elevation data?


texts: Got different cell sizes from rasterised data?

texts:
I'm trying to rasterize the Land cover data of Nigeria which was obtained from FAO database. During the process of rasterizing the land cover data, I used the following information which gave different cell sizes between the input cell size and output size (measured).

From polygon to raster, I employed;

Input feature: selected land cover data layer (nga_gc_adg)
Value field: GRODCODE
Output raster dataset: H:\Arcgis\defaul.gdb\nga_gc_adg_polygonTo raster1
Cell assignment type: CELL CENTRE
Priority field (optional): None
Cell size (optional):0.038.

Why, instead of the same 0.038km used as the cell size, when I measured  the smallest cell size of the rasterized data, do I get much higher cell size (4.24km) compared to the input cell size (0.038km)?


texts: How to retrieve the type of CRS (geographic or projected)?

texts:
Is there a way with GeoTools to know whether a CRS is geographic (in degrees) or projected (in meter) ?

I have tried empirically with the following code (based on a check of the CRS unit) but I am not sure at all if it is going to work for all CRSs... Is there a better and more correct way ?

public enum CRSType { GEOG, CARTO, UNKNOWN }

private static CRSType getCRSType(Unit<?> unit) {
    if(unit == null) return CRSType.UNKNOWN;
    switch (unit.toString()) {
    case "": return CRSType.UNKNOWN;
    case "°": return CRSType.GEOG;
    case "deg": return CRSType.GEOG;
    case "dms": return CRSType.GEOG;
    case "degree": return CRSType.GEOG;
    case "m": return CRSType.CARTO;
    default:
        LOGGER.warn("Unexpected unit of measure for projection: "+unit);
        return CRSType.UNKNOWN;
    }
}

public static CRSType getCRSType(CoordinateReferenceSystem crs) {
    return getCRSType(CRSUtilities.getUnit(crs.getCoordinateSystem()));
}


texts: How to set different map layers for different zoom levels using ArcGIS JavaScript API?

texts:
I have created a layer using "clip" tool in ArcMap.

I have another layer in TOC too. whichis the same area in "clip" Layer with more polygons.

In my web map application (javaScript API) Is there any way to set these layers in a way that when I zoom in I get more details like image2 and by zooming out I get polygons in image1.


texts:
I have a SQL expression that I use in Make Feature Layer in Model Builder ("HOSPITAL" = '% HOSPITALNAME%'). The expression makes a selection in a OD Cost Matrix that a patient is patient only at one hospital and not all hospitals as an OD Cost Matrix would normally do. HOSPITAL is a column in orgins in add location and HOSPITALNAME is a column in the DESTINATION to add location.
What I wonder is how this expression will be in a python script?
I've tried something like this in python in the expression field in arcpy.MakeFeatureLayer "'HOSPITAL' = '" + HOSPITAL NAME + "'". But it does not works. I don't get any error message but it only runs OD Cost Matrix normally.


texts: Executing async method inside relaycommand using ArcGIS Pro SDK and MVVM?

texts:
I'm developing a "search tool" dockpane in the ArcGIS Pro SDK using MVVM. I have a button on my dockpane hooked up to a relaycommand in my viewmodel, which trys to call an async method in my model to return a list of results:

    /// <summary>
    /// Search results
    /// </summary>
    private ObservableCollection<string> _searchResults = new ObservableCollection<string>();
    public IList<string> SearchResults
    {
        get { return _searchResults; }
    }

    /// <summary>
    /// Search button command
    /// </summary>
    private ICommand _searchCommand;
    public ICommand SearchCommand
    {
        get
        {
            _searchCommand = new RelayCommand(
                async() =>
                {
                    SearchResults.Clear();
                    var results = await DockFindModel.SearchAsync(_selectedSearchableLayer, _searchString);
                    foreach (var item in results)
                    {
                        SearchResults.Add(item);
                    }
                    //notify results have changed
                    NotifyPropertyChanged(() => SearchResults);
                },
                () => bAppRunning); //command will only execute if app is running
            return _searchCommand;
        }
    }

However I get the following exception when the relaycommand tries to execute:

  An unhandled exception of type 'System.AggregateException' occurred in
  mscorlib.dll

  Additional information: A Task's exception(s) were not observed either
  by Waiting on the Task or accessing its Exception property. As a
  result, the unobserved exception was rethrown by the finalizer thread.

I've tried a number of things in this thread to try and resolve the issue with no luck. I'm sure running aync methods in a relaycommand is pretty fundamental to the ArcGIS Pro SDK as ESRI recommend you use MVVM, plus a lot of the supplied methods must be ran asynchronously.

Anyone have luck with this?


texts: Converting feature set inputs into feature class or feature layer for Python Script tool

texts:
I want to accept feature set as input and convert it into feature layer or feature class so that I can use the input in select by location in the python script, as the select by location don't accept feature set as input.

For the graphicFeatureSet input I am drawing a polygon and for input_featureSet I am accepting to select featureSet within the drawn graphic layer.
PS: I don't want to use the input type to feature class/layer it should be a feature set.

graphicFeatureSet = arcpy.GetParameterAsText(0)
input_featureSet = arcpy.GetParameterAsText(1)

arcpy.MakeFeatureLayer_management(graphicFeatureSet, "lyr")
arcpy.MakeFeatureLayer_management(input_featureSet, "input_lyr")

arcpy.SelectLayerByLocation("input_lyr","WITHIN", "lyr")


texts:
I have created a script (see below) that works ok but only in case if all parameters are required inputs. Quite simple, Roads(line) and Trackts(polygon) feature classes , first buffer Roads and then clip to Tracts. There are more feature classes in Zion.gdb but I am just doing a test script to figure out optional input parameters.

What I need is that user should be able to select the layer they want to buffer, the layer they want to clip, the output location/name of the clip, and the buffer distance. All of it has to be optional.
The problem is when I set for example Tracts and Tracts_Clip feature classes to optional in Script Parameters and do not apply them in script execution because I want only Roads buffer, the script does not work.

import arcpy

from arcpy import env
    # set environment settings
    env.workspace = "C:/data/Zion.gdb"
    env.overwriteOutput = True

    # Script arguments

    Roads = arcpy.GetParameterAsText(0)
    if Roads == '#' or not Roads:
        Roads = arcpy.GetParameterAsText(0) # provide a default value if unspecified

    Tracts = arcpy.GetParameterAsText(1)
    if Tracts == '#' or not Tracts:
        Tracts = arcpy.GetParameterAsText(1) # provide a default value if unspecified

    Tracts_Clip = arcpy.GetParameterAsText(2)
    if Tracts_Clip == '#' or not Tracts_Clip:
        Tracts_Clip = arcpy.GetParameterAsText(2) # provide a default value if unspecified

    Buffer_Distance = arcpy.GetParameterAsText(3)
    if Buffer_Distance == '#' or not Buffer_Distance:
        Buffer_Distance = arcpy.GetParameterAsText(3) # provide a default value if unspecified

    # Local variables:
    Tracts_Clip = arcpy.GetParameterAsText(2)
    Buffer_Distance = arcpy.GetParameterAsText(3)
    Roads_Buffer = "Roads_Buffer"

    # Process: Buffer
    arcpy.Buffer_analysis(Roads, Roads_Buffer, Buffer_Distance,"FULL", "ROUND", "ALL", "", "")

    # Process: Clip (1)
    arcpy.Clip_analysis(Tracts, Roads_Buffer, Tracts_Clip, "")


texts: Transverse Mercator Projection formula for WGS84 to UTM zone40N (32640)

texts:
I want to convert Lat Long (in WGS84) to UTM zone40N (EPSG::32640) in Excel. So I'm creating a macro for it by following the EPSG Guidance Note 7, part 2.

I'm following the JHS Formula for Transverse Mercator mentioned in page 52 and the examples after it. It seams easy to just calculate the constants for my area and put them in the formula.

The info I have now is:
(reference )

Latitude of Origin (Lat0) = 0
Longitude of Origin (Lng0) = 57
Scale of factor = 0.9996
False Easting= 500,000
False Northing= 0

To calculate the first constant I have the formula

QO = asinh(tan(LatO)) – [e atanh(e sin(LatO))]

And the later constants can be calculated from it, But to calculate Q0, I don't know the values 'a' and 'e'.

In the example 'a' and 'e' are calculated as a parameter of the ellipsoid used for Projected Coordinate Reference System OSGB 1936. I'm assuming the values will be different for my projection system.

Please help me from where can I get these values.

Also let me know if I'm moving in the right direction or not, or if there's another way for achieving my goal. God bless the libraries like ArcPy and GDAL for saving me from calculations like these for so long.


texts: Creating polygon of gaps/holes from several shapefiles using ArcGIS Desktop?

texts:
How can I create polygon of gaps/holes if I have several shapefiles.
Do I really have to merge all layers, it it big amount of data (NEW YORK roadbed, buildings, water,parks etc)

I have ET geowizards to clean polygons or create gaps but I can only choose one layer. Merging all layers takes so long time.

maybe it is another way? erasing, clip with an background layer?.

Union with no gaps allowed worked but i need to be able to symbolize the gap polygon and i need to be able to select individual gap polygons.

any ideas?


texts: ArcGIS Android export tile cache for higher levels (only level 9 in samples)

texts:
I need to cache a few tiles for the offline part of my application. Though in the samples provided for Android the following server link has been used for exporting local tile chache:

http://sampleserver6.arcgisonline.com/arcgis/rest/services/World_Street_Map/MapServer

the sample itself can be found here:

https://developers.arcgis.com/android/sample-code/export-tile-cache/

Where only levels 1-9 are supported.
How can I download for higher levels?
Or, what is the ESRI server link to export higher level tile packages?

P.S. I have already tried other REST services


texts: Relative Path not working in ArcMap and ArcObjects

texts:
I'm using ArcMap 10.0 and Windows 7 64 bit. I'm trying to set a relative paths for my MXD file. My MXD is located at C:\workspace\rastertest.mxd and the data is located at C:\workstation\data\testdted.dt2.

When I set Relative Path to true in ArcMap's Document Properties or in ArcObjects using IMapDocument.SaveAs(true, true) the new raster data source will point to C:\Users\tester\AppData\Local\Temp\data.

Do you know what could be causing this? I would like to eventually move around the MXD file and Rasters and keeping the same directory hierarchy.

Also, when I run ArcMap 10.0 on a Windows XP machine the relative paths work correctly.


texts: Outputting Polygon Neighbors into special format using ArcPy?

texts:
I need Polygon Neighbors to give me data in a specific format. Every polygon should have its own line in the output. Take this, for example(it's obviously not ArcGIS formatted, but it serves the purpose of an example):

This should give the following output, in text or excel:

1 2 5
2 1 3 5
3 2 4
4 3 5 6
5 1 2 4
6 4

It's important that the first column be in ascending numerical order. I have over 300,000 lines of this.

I have the input feature class. I can use ModelBuilder and Python.


texts: Creating connected network of roads from overlapping roads using R?

texts:
I am trying to create a dataset of roads data to do shortest path analysis. To do this, my roads data needs to be connected: it should be possible to jump from one road to another. I am planning to use the R library shp2graph for this. This library contains a function nt.connect that finds the largest connected network. However, this function finds that my data hardly overlaps, despite appearing to overlap on visual inspection. The function returns the following image to visualise these separate networks:

How can I connect these overlapping pieces of road to each other?

What I've tried:

Buffering the roads network. This returns a SpatialPolygon object which nt.connect does not accept: it wants a SpatialLines object. If I would convert the SpatialPolygon object to a SpatialLines object it would take the boundary of the Polygon, which would result in two roads right next to each other. I am worried this will mess with the subsequent analysis.

Some notes:

I am using OpenStreetMap data. The dataset (a subset of my full dataset) can be downloaded here: https://www.dropbox.com/sh/qu2fg3f9fy5tk7k/AABjv81Vv1R9c0RHuMJwIHH2a?dl=0
This is only a (very small) portion of the full roads network I want to use. Manually making sure they overlap is therefore not feasible
This library I am planning to use uses R. My work needs to be fully reproducible and this is difficult when working with a GUI.
Some pieces of the road are true "orphans": they do not overlap with the larger network. I would prefer if these can be incorporated into the main network, for example by drawing the shortest possible line to the full network. I would also be fine with just dropping these orphans.

Below is the code to check that the network is connected:

require(sf)
require(sp)
require(shp2graph)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
roads <- read_sf("../GIS Stack example", "roads")
roads_sp <- as(roads, "Spatial")
roads_connected <- nt.connect(roads_sp)


texts: How do you manage existing GIS data and allow for integration with an Asset management system?

texts:
I am using ESRI based GIS software, Postgresql/PostGIS/ArcSDE DB and we have a Mincom Ellipse asset management system.

Currently our all our spatial information regarding asset sites have been recorded as points, this has been versatile for the fact that it accommodates mapping at various scales. Now that we are integrating our Asset Management System with our GIS database the asset management guys want the GIS features to reflect the structure e.g A building footprint as a polygon instead of a point.

My question in terms of the spatial data management. Should I be maintaining two sets of data? One for the asset representations and then one for various mapping tasks?

Thanks
DB


texts: Overlapping symbol in print layout in QGIS

texts:
I am working on water network map in QGIS. I've a layer with different types closely spaced of symbols.

For exporting the map, when i zoom out the map some symbols overlapped with each other & the map looked messy.

Is there any way to show all symbols clearly in the map?


texts: Clutter data - where to find for BC and how to convert to .grc

texts:
I would like to find a source of 30m or 90m clutter data for BC, Canada. Can anyone point me in the right direction? Preferably a free source.

Is landcover synonymous with clutter? when downloading clutter data, how do I differentiate between grid elevation data and clutter data?  Is .lcv clutter data? could clutter data be in the form of a geotiff?  A basic primer on this would be great if anyone has any links.

I am using a radio propagation planning tool and it accepts clutter in a number of formats including .grc. I can use a tool called Global Mapper to convert to .grc if needed. Are there other/better tools to convert clutter data?


texts: Bulk Import rasters into GRASS QGIS plugin

texts:
I'm trying to batch import files into the GRASS (v. 7.4.0) plugin in QGIS 2.18.18. I'm sure this is pretty basic, but I'm not sure quite how to do this. This exact question was asked before here, but it was for an earlier version of QGIS. Does anybody know how to to do this for 2.18.18?

EDIT: The rasters are different folders, all contained in the same parent folder (with all the same name no less...). I could add all of them in QGIS by dragging groups of folders into the window, is there anything like this for GRASS?


texts: Find nearest distinct features to a point

texts:
I am trying to select a single distinct features closest to an XY point in PostGIS.

So I can find a single record for each distinct code by using:

Select Distinct ON (code) * FROM table;

Which returns 100 records, out of 1mill, one for each distinct code.

I can then use the PostGIS <-> operator to find the closest feature:

SELECT name, gid
FROM geonames
ORDER BY geom <-> st_setsrid(st_makepoint(-90,40),4326)
LIMIT 1;

So I need to use both of these queries together to find a feature of each distinct code but the closest one to the point.

How can I use both queries together?

I have tried:

SELECT 1 FROM table WHERE EXISTS
 (SELECT DISTINCT ON (code) * FROM table) ORDER BY geom <-> st_setsrid(st_makepoint(-90,40),4326)
LIMIT 1;

Also this fails:

SELECT * FROM table WHERE code IN
(SELECT DISTINCT ON (code) * FROM table) ORDER BY geom <-> st_setsrid(st_makepoint(-90,40),4326)
LIMIT 1;

None of which works.


texts: How to gain the facility of snapping when multiple points are selected

texts:
I have created a fishnet label(square 961 points) centering a point in other layer in arc map but the fishnet is deviated from the intended center.
Now i want to place the midpoint of labels on the intended center by selecting all the labels.
But when i select all label points, i do not get the facility of snapping whereas i put the midpoint of the labels on top of the intended center accurately.
Whenever i zoom the map i found midpoint of fishnet labels and intended center did no match.

                                Present situation

                                Wanted situation

Both are in the same coordinate system.
How can i get rid of this situation .

Thanks in advance


texts: Filter features with multi-table joins in ESRI Silverlight/REST API

texts: We have a rather involved data warehouse (SDE 9.3.1 / Oracle 11g) that contains a single related point feature class. Our (analytical) application should allow users to construct complex, ad hoc queries against the data. In every case there will be joins across numerous (but not always the same) tables; only matching features should show on the map. The REST & Silverlight APIs support a simple where clause, which seems to limit our options. Can you recommend a strategy for surmounting this limitation?


texts: How applications like Global Mapper, ArcMap and ArcScene support showing high resolution geotiff images?

texts:
Considering my question in Game Developement:

How to create textures of big images in opengl (bigger than the MAX_TEXTURE_SIZE)

it seems that, there's no way to show an image bigger than the MAX_TEXTURE_SIZE in OpenGL (I'm not sure yet, looking forward to more answers in Game Developement)

In my case MAX_TEXTURE_SIZE = 8192 but my image is a 16997x15931 which can be shown in openGL based on my current knowledge about this API.

Just tobe curious, and maybe they're using a toolkit or API which is not free

What way programs like Global Mapper, ArcMap and ArcScene manage to show high resolution geotiff images for example bigger than 15000x15000 with high performance?

What toolkit, API, 3rd party library do they use for display tasks and 2D or 3D rendering in their programs?

I highly doubt that image-pyramids is the concept related to my question?
if it's true beyond answering my question about APIs please introduce me a link in GDAL documentation that has an excellent tutorial about pyramids?


texts: Selecting MAX values within records with the same name

texts:
I'm working in ArcMap 9.3 and I've never dealt with that kind of queries...

I need to Select the MAX value from column [A] for each set of records in the table which have the same value in column [B]

For example

[A]  [B]

100  aaa

101  aaa

102  aaa

999  aaa

100  bbb

101  bbb

102  bbb

999  bbb

so the selection would be record#4 (where [A]=9999 which is the biggest value of all 'aaa') and record#8 (where [A] = 999 which is the biggest of all 'bbb')

Is it possible to do that using queries and subqueries when working with databases?

Any help appreciated!
Dan


texts:
I try to insert a line between two lines. End points of lines are in intersection of lines. I use positionAlongLine in arcpy but it interestingly gets coordinates from the first points of referenced lines even if I use shape.positionAlongLine(0.05).lastPoint Also I get 5cm distance but it generated very long lines. The lines created like the picture. Red lines are generated, blue lines are referenced. The result is oposite of what I want. How to acquire a line using positionAlongLine?

import arcpy
arcpy.env.workspace = "C:\Users\Asus\Documents\workspace\PROJECTS\output\output_data.gdb"
coords = []
point = arcpy.Point()
array = arcpy.Array()
featureList = []
cursor = arcpy.InsertCursor(r"anahtar_line")
feat = cursor.newRow()

with arcpy.da.SearchCursor("line", ["OBJECTID", "SHAPE@"]) as scursor:
    for row in scursor:
        shape = row[1]
        end = shape.lastPoint
        point.X = end.X
        point.Y = end.Y
        array.add(point)
        d = shape.positionAlongLine(0.05).lastPoint
        point.X = d.X
        point.Y = d.Y
        array.add(point)
        polyline = arcpy.Polyline(array)
        array.removeAll()
        featureList.append(polyline)
        feat.shape = polyline
        cursor.insertRow(feat)
del feat
del cursor


texts: Adding a WKT specified coordinate system to ArcGIS?

texts: I've got a WKT specified coordinate system, how do I add this to ArcGIS?


texts:
I'm doing a temporal analysis of NDVI and temperature over 17 years, necessitating data from Landsat 5,7, & 8. I'm interested in using level-2 Landsat data available from https://espa.cr.usgs.gov/ordering/new/. I would be able to directly download NDVI as well as brightness temperature, which would save me a lot of time. However, the USGS is explicit that the algorithm used to calculate surface reflectance for Landsat 8 is not the same as the algorithm used for previous satellites.

https://landsat.usgs.gov/landsat-surface-reflectance-high-level-data-products

My question is: are the level-2 products from Landsat-8 comparable on an apples-to-apples basis with Landsat 5 or 7 level-2 products, or does this kind of comparison require some sort of adjustment to the data? I haven't been able to find an answer to this question online.


texts: "ImportError: No module named fiona" (post-installation of package)

texts:
I'm in need of the Fiona and Shapely modules to develop a QGIS 2 plugin, and thus downloaded them from http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal and installed them. I have no trouble importing shapely, but importing fiona gives me "ImportError: No module named fiona"

I am wondering if there's a step I'm missing besides the installation of the correct version (for 32-bit Windows and Python 2.7, using the standalone installer for QGIS) before being able to import them via the QGIS Python console.


texts:
I have a shapefile and their X and Y coordinates' values are in Decimal Degrees. I want it to convert it in meters.

I also tried to use Calculate Geometry tool. But as you might see it from below screenshot, I am only getting a single unit, i.e., Decimal Degrees under the 'Units' drop down.


texts: Make a condition in QGIS with other layers

texts: I want a QGIS composer form to change its color according to whether the value of an attribute of a layer matches a value, but for this I want to choose the attribute of a layer that is not the current.


texts: Knowing if error occurred while loading WMS layer in OpenLayers?

texts:
I already asked this question in Stack Overflow. I don't want to spam by asking again here but I think the question fits better in GIS SE. I am sorry for the replication. Feel free to delete one of the two questions.

Let's start!

I am working with OpenLayers 2. I am using the events "loadstart" and "loadend" from here.

When I use OpenLayers.Layers.WMS, i want to have some kind of callback when there is a loading error. For example, when I enter an incorrect SLD filter, I want to have some feedback about the layer not being loaded correctly. Unfortunately, "loadend" doesn't provide this kind of error detection.

I already found this question where someone else has the same problem.

Is there even an event that signals this kind of error?

Is there something like a "layerfail" or "layererror" event?

How could I implement this?


texts: `pyproj` initialization error in PyCharm but not in cmd?

texts:
I am currently struggling to get pyproj to initialize when running in PyCharm. Oddly enough, I am able to get it to work in cmd.

The version of pyproj I am using is 1.9.6, and the version of proj4 I am using is 5.2.0.

The following runs without error using the Python interpreter:

import pyproj
pyproj.Proj({'init': 'epsg:26910'})

However, when running from PyCharm the following error is produced:

Traceback (most recent call last):
  File "C:/Users/frankbr/.PyCharmCE2019.2/config/scratches/scratch.py", line 5, in <module>
    pyproj.Proj({'init': 'epsg:26910'})
  File "C:\Users\frankbr\.conda\envs\pyfor_env\lib\site-packages\pyproj\__init__.py", line 362, in __new__
    return _proj.Proj.__new__(self, projstring)
  File "_proj.pyx", line 129, in _proj.Proj.__cinit__
RuntimeError: b'no arguments in initialization list'

Due to other dependencies, and the fact that it runs in cmd (and Travis tests) without error, I am hoping to find a solution that only requires modification to PyCharm.


texts: How do I get shapely polygon intersection area?

texts:
In an example case, at debug point, y=40000, x=1500, z=0. Small polygon is totally inside big and I thought z should be 1500. contains doesn't find the case at all.

for big_polygon in big_polygons:
    for small_polygon in small_polygons:
        if big_polygon.intersects(small_polygon):
            y=big_polygon.area
            x=small_polygon.area
            z=big_polygon.intersection(small_polygon).area
            #debug point


texts: This page [https://cartodb.com/cartodb-for/twitter-maps/] seems to indicate a live tweet map is possible with CartoDB, but the interface forces me to define a time span and doesn't allow me to define beyond the current date. Does anyone know how to create a Tweet map in CartoDB that will update live as new tweets meet the search criteria?


texts: Is there anyway in FME to create spatial features from data in Microsoft CRM?

texts:
Is there anyway in FME to create spatial features from data in Microsoft CRM Online?

We have various CRM entities, which we'd like to put on the map. Some have addresses, some have UK OS NGR, some even have lat/long values.

I've checked the list of FME formats but can't find anything relating to Microsoft CRM

Somebody suggested going directly to the SQL Server tables, but this isn't possible as we're using CRM Online, so database connection


texts:
I'm having trouble trying to chart multiple regions onto one histogram.
var urban1 = urban.first().geometry();
var pine_forest1 = pine_forest.first().geometry();

var regions = ee.FeatureCollection([
  ee.Feature(urban1, {label: 'urban'}),
  ee.Feature(pine_forest1, {label: 'pine_forest'}),
]);

var hist = ui.Chart.image.histogram({
  image: may_filtered_clip,
  region: regions,
  scale: 20
}).setOptions( {
  series : {
      0 : {color : 'red'},
      1 : {color : 'green'}
    }
  });

print(hist);

This is the result of the code:

And this is basically what I'm aiming for:


texts: How to calculate statistics for several sub-catchments (polygons) in the same time?

texts: I have a shapefile of groundwater management zone which is composed of 40 sub-catchments (40 polygons). I'd like to calculate the max., min. and mean elevation for each sub-catchment. Would you please let me know if there is any way to calculate all of these parameters (for all of the sub-catchments) using the shapefile and the DEM?


texts:
I have a database of latitude, longitude and radius in the following form:

48180759,11518950,19.0

Both latitude and longitude are expressed in degrees but then multiplied by 1e6 (or 1000000). Radius, on the other hand, are expressed in meters. What I am trying to do is convert lat/lon into meters so that I can use geometric functions more accurately, especially the one that looks for intersects (overlaps between points).

EDIT:
Unfortunately, I wasn't able to find a complete answer. At least, not in this web site. I will use Python for this task. As for the precision, I would like to have as best as possible (a few centimeters deviation at most).

Any ideas?


texts:
I'm trying to implement creating custom geographic object by GeoDjango. In model I have this field.

mpoly = models.MultiPolygonField()

Then I create form

 from django.contrib.gis import forms as geoforms
 from world.models import Country

 class CountryForm(geoforms.ModelForm):
    class Meta:
        model = Country
        widgets = {
            'mpoly': geoforms.OSMWidget(attrs={'map_width': 800, 'map_height': 500}),
        }

and rendered it by

<head>{{ form.media }}</head>
<body>
   <form>
      {{form}}
      <button class="btn btn-primary" type="submit">Save</button>
   </form>
</body>

Full page code is here

After drawing polygon on map I get error in MapWidget.searializeFeatures function.

Uncaught TypeError: a.Y is not a function
at bp (https://cdnjs.cloudflare.com/ajax/libs/ol3/3.20.1/ol.js:459:803)
at Zo.k.Te (https://cdnjs.cloudflare.com/ajax/libs/ol3/3.20.1/ol.js:464:27)
at Zo.k.od (https://cdnjs.cloudflare.com/ajax/libs/ol3/3.20.1/ol.js:431:47)
at MapWidget.serializeFeatures (http://localhost:8000/gstatic/gis/js/OLMapWidget.js:227:69)
at D.<anonymous> (http://localhost:8000/gstatic/gis/js/OLMapWidget.js:89:22)
at D.b (https://openlayers.org/en/v4.0.1/build/ol.js:41:394)
at Dc.b (https://openlayers.org/en/v4.0.1/build/ol.js:44:213)
at D.k.Be (https://openlayers.org/en/v4.0.1/build/ol.js:47:617)
at D.k.push (https://openlayers.org/en/v4.0.1/build/ol.js:48:41)
at Qt.k.Nd (https://openlayers.org/en/v4.0.1/build/ol.js:634:301)

I can't find anything meaningful about this error.

What is the reason for this error, and how can I fix it?


texts: Using QGIS to plot locations on map of US by zip code?

texts:
I am new to QGIS.  I just installed QGIS 2.2 on my MAC Airbook with Mavericks.

How do I simply plot a few locations on a map of the US?

All of the tutorials I can find on the web are for older versions.  Version 2.2 seems to be entirely different in terms of look and feel.  For example, one tutorial tells me to load the plugin "edit any layer".  Under 2.2 no such plugin seems to exist. I simply need to know how to get map loaded and then how to plot some excel or CSV data.


texts: ArcMap Add-In Not Updating To Latest Version

texts:
I'm having an issue that I haven't experienced before with ArcMap Add-Ins and I'm hoping someone out there has experienced this and knows of a solution.

We have an add-in build in VB.NET that is used by multiple users in my office. It is stored in a central location that we all have access to and each person who needs to use it has that folder added to their Add-In Manager. That way any needed updates can be made and everyone gets the latest copy every time they open ArcMap.

That has worked well until today. One user in particular can't seem to get the updated version. Even though the Add-In Manager shows the correct version when he goes to use the add-in it is clearly an old version because it doesn't have the latest change made to the code. Everyone else (pointed to the exact same folder) gets the latest version. We even removed the folder from the Add-In Manager and manually installed the add-in. Same thing. He always gets an old version no matter what we do.

Is there some type of cache I'm not aware of that might be holding onto the old one and tricking ArcMap into thinking it has all the updates?


texts: How can I make the UTMs of different datasets uniform?

texts:
I have wind datasets for every U.S. state and want to merge them. I found out that these datasets use different UTMs. So I want to make the UTMs uniform. How can I do that? This can not simply be done by a projection, can it?

The metadata says:

  This shapefile was generated from a raster dataset with a 200 m
  resolution, in a UTM zone 13/14/15/16 etc., datum WGS 84 projection
  system.

This is my data source: http://www.nrel.gov/gis/data_wind.html

I am using ArcGIS 10.2.2


texts: Only show features, related to a selected feature in another layer (QGIS)

texts:
I have: ways (tracks) and points (way signs). The tracks are outgoing from the way signs.

This should be the visible result:
The ways with name 'A' are visible (coloured) because the related way sign 'A' is selected. This is the result I'm looking for:

The question is: How do I show or high light only the tracks, related to the selected way sign?

What I did:
With this rule in styling settings I can realized the wanted result with one layer (waysigs OR tracks): I select the way sign /tracks and it becomes red, because I defined it as styling rule.

Then I tried this for the tracks (without understanding, what I'm doing) related to the example-syntax in expression builder Help section¹:

    is_selected(get_feature('waysigns', 'fid', "name"), 'waysigns')

...but without the expected result...

  ¹
  s_selected(get_feature('streets', 'name', "street_name"), 'streets')
  → True if the current building's street is selected.

Why do I need this: In my project I have thousands of tracks and some way signs. A lot of the tracks are crossing each other and it is impossible, to see, which track comes from which way sign. To have a better overview, which tracks are related to which way sign, I simply would like to select the way sign, and get its related tracks.


texts: problem with feeding input files into gdal_merge.py

texts:
While merging the files using gdal_merge.py the following code works:

subprocess.call([sys.executable,gmerge,'-o','C:\\r.tif','-of','GTiff','D:\\a.tif','D:\\b.tif'],shell=True)

However, when the input files are numerous, each file can not be inserted separately. In this case, the following code does not work:

subprocess.call([sys.executable,gmerge,'-o','C:\\r.tif','-of','GTiff','D:\\*.tif'],shell=True)

I did not know how to input the input files.
Any idea is welcomed.


texts: Retrieve geometry type (point/line/polygon) from OL3 feature

texts: How do I tell whether an OL3 feature is a point, line, or polygon? I need to highlight features on mouseover, and the style required is different for points, so I need to be able to tell which geometry type a feature is, so that I can apply the correct style.


texts: Imported GPX files into QGIS ver 2.0.1-fields missing?

texts:
Just imported GPX file into QGIS.Opened attribute table.
where are the following fields? Latitude,Longitude,,symbol,and time . Am I missing something??I used a Garmin Oregon 300 to collect the data. Verified that all the listed fields have info in them using other both Mac and Windows GPS programs.


texts: Calculating Sinuosity Using ArcGIS Field Calculator

texts:
Trying to calculate sinuosity via Field Calculator (ArcGIS v 10.2.2); VBA

Pre-Logic Script Code:

Dim dblLe as double
Dim dblLf as double
Dim dblS as double
Dim pCurve as Icurve
Dim pFPoint as IPoint
Dim pTPoint as IPoint
Set pCurve = [Shape]
Set pFPoint = pCurve.FromPoint
Set pTPoint = pCurve.ToPoint
dblLe = Sqr((pFPoint.X - pTPoint.X)^2 + (pFPoint.Y - pTPoint.Y)^2)
dblLf  = pCurve.Length
dblS = dblLe/dblLf

Variable =

dblS

I know ESRI has a python script put I am try to run via Field Calculator.  Anyone see my error?


texts: Save edits to polygon without altering existing

texts: I am trying to edit a single polygon and save it as a new file without altering the original polygon in QGIS 3.4.


texts:
How to raster clip in ArcMap? But I need 1 raster in outside area. Because 2 raster map.

Example:

Ras1=Temperature.tif
Ras2=Error map.tif
Clip command input raster=Ras1
             Output Extent= Ras2

Just I need clip outside area=reserve clip command? I can't use Con and Setnull, IsNull.


texts: Problem changing CRS to enable clip function

texts:
I have two layers, one point file and one polygon.  The point file is in WGS 84 and polygon in NAD83-UTM10.  The project CRS is NAD83-UTM10.  They project fine with on-the-fly projection turned on.  I want to clip the point file with the polygon, and understand that all layers need to be in the same CRS to enable clipping or other ftools functions. When I change the CRS of either layer to the other (changing the point file to NAD83 for example), it projects far away from where it was (ie. doesn't line up).  This is the same if I right click to change the CRS or if I "save as."  Same no matter which CRS I select.  I'm running QGIS 2.0.1 Dufour on Mac.

Any help would be greatly appreciated!  Thanks!


texts: How to operate maximum value compose of some raster image?

texts:
I want to find the maximum/minimum value from a serials images at each pixel location. These images have the same resolution and projection.
such as  Z（m）= max(Z1,Z2,Z3....) where Z indicates a raster file.

I know ArcMap raster calculator tool contains a Function 'con'. Are there any rapid methods in Python or GDAL?


texts: Hosting a few shapefiles so multiple users could use and edit the same map

texts:
Is there a way to synchronize a QGIS map with vector layers so that 3 or 4 users can use and edit the same file?   I have a plug-in that updates shapefiles on a layer and it would be great to synchronize this on a cloud service.   I have tried QGIS Cloud and it doesn't seem to work well with vector layers.

http://support.qgiscloud.com/kb/index.php


texts:
I know this has been asked and answered here but since I don't have enough rep to ask a clarifying question to the person who answered it I am having to post a completely different question.

I am trying to perform zonal stats on a geotiff with 6 bands (RGB and three thermal and NIR bands). Like the post below it only calculates the stats on the first band.

How does ArcGIS calculate zonal statistics with multi-band rasters?

In the screen shot in the above post, it appears that you can select a specific band when you select the "Input value raster", however mine doesn't look like this. I simply get the Raster Datasets rather than the Raster Bands.

So my question is, how do I get the individual bands so I can calculate zonal stats on each band in this way? Is it that I have to somehow export single bands as new rasters - if so how do you go about doing this?


texts:
I need the 'concave hull' or alphashape of a set of points. Also, I'd like to have it with holes/donuts/islands in the middle. See my attached image. The outer hull is being calculated fine. But the large islands in the middle are generated as part of the hull. I want to get rid of that.

My code (partially grabbed from tutorials all around) can be found here: https://gist.github.com/willemvanopstal/f199c2f839b4de2abadfafde450c45ed

I expect it to be a simple mistake from my side.. Maybe in generating the resulting polygon from the isolated edges; (edge_points are the edges which have passed the alpha-test)

m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points

Is there anybody who knows a fix for my issue?


texts: Esri uniform meters grid to geographical coordinates

texts:
I have a bathymetry for which I have a grid file:

Projection    MERCATOR
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters
 100  0  0.0 /* longitude of central meridian
 -41  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)

and a bathymetry file itself:

ncols 12116
nrows 15391
xllcorner 4795705.6
yllcorner -5915585.8
cellsize 250
NODATA_value -9999

I understand that 250 m grid is valid only for -41 latitude of true scale, so I need to scale for other latitudes. However, I don't understand what means xllcorner/yllcorner in this case. How can I convert them to latitude and longitude to use as a reference point? Longitude of central meridian (100deg) doesn't seem to be what I need, because the bathymetry is for New Zealand.


texts: Changing coordinate system in which WMS is served using ArcPy?

texts:
I am trying to change the coordinate system of a WMS.
what i have tried is:

prjFile = lokaal + "opmeting_straatkolk.prj"
sr = arcpy.SpatialReference(prjFile)
mxd2 = arcpy.mapping.MapDocument("CURRENT")
for DF in arcpy.mapping.ListLayers(mxd2):
    DF.spatialReference = sr
    arcpy.RefreshTOC()
    arcpy.RefreshActiveView()


texts:
I am facing troubles geocoding the CSV file in QGIS 3 with MMQGIS (OpenStreetMap).

Images below are showing my workflow and the problem that I am encountering.

I am getting the following error:

Failure connecting to nominatim: <urlopen error [SSL:
CERTIFICATE_VERIFY_FAILED] certificate verify failed (_ssl.c:852)>: 0
addresses geocoded

I'm using a Mac with OS 10.14.1 (18B75) and my Internet connection is stable.

These are my QGIS specifications:

Any ideas?


texts: Export UTM coordinates with sub-meter precision from a shapefile in ArcGIS Pro?

texts:
I have a series of points in ArcGIS Pro that were collected with sub-meter precision. I need to export UTM coordinates for each point so that I can use them for some processing outside of Arc, and I need them to be accurate to 2 decimal places to be consistent with my other data.  I tried using the Calculate Geometry Attributes tool, but the coordinates it outputs are rounded up to meters (and yes, I have set the display in the attribute table to show the correct number of decimal places).

Is it that Calculate Geometry doesn't calculate outputs to that level of precision?  If so, what tool should I use instead (I can do some basic Python, if necessary).

Or... is there something wrong with my underlying data?  I collected the data with a Trimble GeoExplorer handheld (can't remember the exact model), did a differential correction with Pathfinder office, brought the output shapefile directly into Arc, and then projected and transformed it from WGS84 Lat/Long to NAD83 UTM 10N. Is there a step in there that could have truncated my coordinate values?  I'm new to working with data at this level of precision.


texts: Converting EPSG:2284 to EPSG:4326 with pyproj

texts:
I have some points in EPSG:2284 and I'm trying to convert them to EPSG:4326 using pyproj, however I get weird results. When I go to the spatial reference website
I see that there is an endpoint that does the conversion. When I use that instead I get the correct result.

Here is some code illustrating the different results:

from pyproj import Proj, transform
import requests

if __name__ == '__main__':

    # points I want to convert from EPSG:2284 -> EPSG:4326
    targets = [
        (3669486.63386,11293013.10592,709.08004),
        (3669559.13811,11292972.72831,711.60055),
        (3669639.51308,11292851.48264,712.22258),
        (3669800.62304,11292949.38118,714.75766),
    ]

    for northing, easting, up in targets:
        print "Converting", easting, northing

        # Transform using pyproj (gives wrong answer)
        WGS84 = Proj(init='EPSG:4326')
        inp = Proj(init='EPSG:2284')
        x, y = transform(inp, WGS84, easting, northing)
        print x, y, "(with proj4)"

        # Transform using spatialreference.org URL (gives correct answer)
        url = "http://spatialreference.org/projection/?json=%7B%22type%22%3A%22Feature%22%2C%20%22geometry%22%3A%7B%22type%22%3A%22Point%22%2C%20%22coordinates%22%3A%5B{lon}%2C%20{lat}%5D%7D%2C%22properties%22%3A%7B%7D%7D&inref=EPSG:{inp}&outref=epsg:{outp}"
        url = url.format(lon=easting, lat=northing, inp='2284', outp='4326')
        x, y = requests.get(url).json()['coordinates']
        print x, y, '(using spatialreference.org URL)'

Does anyone know why they are giving different result and what is wrong with the Pyproj code?


texts:
After following step by step the instructions here to install GDAL: http://sandbox.idre.ucla.edu/sandbox/tutorials/installing-gdal-for-windows

When I run from osgeo import gdal, I get the following:

ImportError: DLL load failed: The specified module could not be found.

How could I fix this?

I use python 2.7.10 :: Anaconda 2.3.0 (64b), Windows 10 (64b)


texts: Saving Processing algorithm output to string instead of file.geojson

texts:
Is there any way to save processing.run("native:algorithm") output directly to string instead of file or temporary layer?

I'm aware of saving output to layer and exporting features using QgsJsonExporter(), however it adds ID to every feature therefore it's not what I'm looking for.

Now I'm going for approach where extra ID is generated for every feature in geojson result is generated. However then I change OUTPUT to file.geojson there are no extra IDs.

vl = QgsVectorLayer(sys.argv[1],"mygeojson","ogr")

def runProcessingNativeBuffer(inputGeoJson, distance):
    processing.algorithmHelp("native:buffer")
    return processing.run("native:buffer",
                         {'INPUT': inputGeoJson,
                          'DISTANCE': distance, 'SEGMENTS': 5,
                          'END_CAP_STYLE': 0, 'JOIN_STYLE': 0,
                          'MITER_LIMIT': 2, 'DISSOLVE': False,
                          'OUTPUT': "TEMPORARY_OUTPUT"})

ats = runProcessingNativeBuffer(vl, sys.argv[2])

vl1=ats['OUTPUT']

features = vl1.getFeatures()
exporter = QgsJsonExporter()
print("GEOJSON",exporter.exportFeatures(features))


texts: Get errors in GeoPy with Python when I load large size data

texts:
I have data size 1.87GB contain 14962000 rows and 9 columns. And I want to apply geopy but I got the errors.
my code:
import pandas as pd
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter

df = pd.read_csv("mydata_name.csv")
df["Geo"] = df["latitude"].astype(str)+ ',' + df["longitude"].astype(str)
geolocator = Nominatim(user_agent="MyApp", timeout=10)

rgeocode = RateLimiter(geolocator.reverse, min_delay_seconds=0.001)

df['address'] = df['Geo'].apply(rgeocode)

When execute my code I get the following errors. I want to get address for all rows. But the code take 5 hours to executing and I get this
RateLimiter caught an error, retrying (0/2 tries). Called with (*('24.807412,46.803666',), **{}).
Traceback (most recent call last):
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\geocoders\base.py", line 360, in _call_geocoder
    result = self.adapter.get_json(url, timeout=timeout, headers=req_headers)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\adapters.py", line 373, in get_json
    resp = self._request(url, timeout=timeout, headers=headers)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\adapters.py", line 401, in _request
    raise AdapterHTTPError(
geopy.adapters.AdapterHTTPError: Non-successful status code 502

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\extra\rate_limiter.py", line 136, in _retries_gen
    yield i  # Run the function.
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\extra\rate_limiter.py", line 274, in __call__
    res = self.func(*args, **kwargs)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\geocoders\nominatim.py", line 359, in reverse
    return self._call_geocoder(url, callback, timeout=timeout)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\geocoders\base.py", line 378, in _call_geocoder
    self._adapter_error_handler(error)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\geocoders\base.py", line 392, in _adapter_error_handler
    raise exc_cls(str(error)) from error
geopy.exc.GeocoderServiceError: Non-successful status code 502
RateLimiter caught an error, retrying (0/2 tries). Called with (*('24.650473,46.669335',), **{}).
Traceback (most recent call last):
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\geocoders\base.py", line 360, in _call_geocoder
    result = self.adapter.get_json(url, timeout=timeout, headers=req_headers)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\adapters.py", line 373, in get_json
    resp = self._request(url, timeout=timeout, headers=headers)
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\adapters.py", line 401, in _request
    raise AdapterHTTPError(
geopy.adapters.AdapterHTTPError: Non-successful status code 502

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "C:\Users\ibrahim\Anaconda3\lib\site-packages\geopy\extra\rate_limiter.py", line 136, in _retries_gen
    yield i  # Run the function.


texts: How to list all names of tile inside a given polygon?

texts:
I'd like to generate a specific set of slippy map tiles to cover certain a certain land mass.

Given a polygon shapefile, is there something that will quickly generate a listing of the z/x/y tile names that sit inside that polygon at a specific zoom level? Ideally, it would allow for buffer around the polygon to give some space.


texts:
am calling a google geocoding api via java, and getting wrong answer ! for example if i run the query below from browser i get the location but if i run it using a rest client i get zero result !!!
is there anything to solve this issue ?

https://maps.googleapis.com/maps/api/geocode/json?sensor=false&address=18+RUE+DE+NIEDECK++++39255+OTTERSWILLER&region=FR&key=AI.....

this is the answer from Google Chrome or Mozilla :

{
   "results" : [
      {
         "address_components" : [
            {
               "long_name" : "18",
               "short_name" : "18",
               "types" : [ "street_number" ]
            },
            {
               "long_name" : "Rue du Nideck",
               "short_name" : "Rue du Nideck",
               "types" : [ "route" ]
            },
            {
               "long_name" : "Otterswiller",
               "short_name" : "Otterswiller",
               "types" : [ "locality", "political" ]
            },
            {
               "long_name" : "Bas-Rhin",
               "short_name" : "Bas-Rhin",
               "types" : [ "administrative_area_level_2", "political" ]
            },
            {
               "long_name" : "Grand Est",
               "short_name" : "Grand Est",
               "types" : [ "administrative_area_level_1", "political" ]
            },
            {
               "long_name" : "France",
               "short_name" : "FR",
               "types" : [ "country", "political" ]
            },
            {
               "long_name" : "67700",
               "short_name" : "67700",
               "types" : [ "postal_code" ]
            }
         ],
         "formatted_address" : "18 Rue du Nideck, 67700 Otterswiller, France",
         "geometry" : {
            "location" : {
               "lat" : 48.7237784,
               "lng" : 7.363452799999998
            },
            "location_type" : "ROOFTOP",
            "viewport" : {
               "northeast" : {
                  "lat" : 48.7251273802915,
                  "lng" : 7.364801780291501
               },
               "southwest" : {
                  "lat" : 48.7224294197085,
                  "lng" : 7.362103819708497
               }
            }
         },
         "partial_match" : true,
         "place_id" : "ChIJXwysfE2nlkcRjwlaKN_gmrE",
         "types" : [ "street_address" ]
      }
   ],
   "status" : "OK"
}

and this is the answer for the same address when i use a rest client :

{
   "results" : [],
   "status" : "ZERO_RESULTS"
}

This is my code , and i am using apache http client version 4.5

    <dependency>
        <groupId>org.apache.httpcomponents</groupId>
        <artifactId>httpclient</artifactId>
        <version>4.5.2</version>
    </dependency>

as for the code here it is :

        String f = "https://maps.googleapis.com/maps/api/geocode/json?sensor=false&address=18+RUE+DE+NIEDECK++++39255+OTTERSWILLER&region=FR&key=googlekey";

    // executing the request
    HttpClient client = HttpClientBuilder.create().build();
    HttpGet request = new HttpGet(f);
    HttpResponse response = client.execute(request);
    String responseBody = EntityUtils.toString(response.getEntity());
    System.out.println(responseBody);


texts: Calculating radius of curvature and banking angle when set of lat, long and elevation are given?

texts:
I have set of lat,long and elevation values with me. Can you help me find the Radius of curve and banking angle.

Formulae may be provided as an answer.


texts: Finding points along lines, given lengths along those lines, using QGIS?

texts:
I have one layer with a network of lines and a second layer with its associated line and its distance along that line.

Is there a way to generate points for all these rows along their associated lines?


texts: Using Dijkstra's algorithm for multiple source to single destination shortest paths

texts:
I'm going over a lecture recording, in it my professor mentions using Dijkstra's algorithm (or a modified version of it) to find multiple-source to single source shortest paths, e.g. finding the closest hospital out of three hospitals to an accident site.

I have no visuals to see what he is doing on the blackboard, but he describes creating a new 'virtual node' that has a node degree equal to the number of sources, with each edge linking to each intended source, and the cost of this node's edges are 0. I guess the process is that the new virtual node becomes the source, and as each of its edges connect to a real sources for zero cost, it would be computing shortest paths from each source in the most efficient manner.

Is this correct?


texts:
I have to get to know the approximate amount of volume of rocks in one area. For that, I need the elevation information to see where slopes or inclinations are located.
I downloaded from the USGS Explorer the GMTED (Global Multi-resolution Terrain Elevation Data) and created a shapefile with the contour lines, where each line represent a particular elevation (I used this tutorial to do so: http://www.qgistutorials.com/en/docs/working_with_terrain.html).

I then sent this information as a .dxf file to a person who wanted to calculate the volume with a program (he told me that it should be a dxf data). But when he opened the data, the elevation lines were just flat and the z values were missing. Does someone know what I have to do now? Is it even possible to get the z data from a DEM the way I did it? Where could I find it? I am just learning QGIS, so everything is pretty new to me. I use the version 2.12.3. –Lyon.


texts:
I'm trying to make an atlas with qgis with a map generating 2 pages (north and south) N and S corresponding to the coverage layer where a field contains the values "N", "S". Every page is formed by the map zoomed on N or S (this is not a problem) and a part where I want the 2 value "N" e "S" changing the font-color: N red and S grey if the atlas zoom the north, N grey and S red if the atlas zoom the south.

I think I can use a conditional expression with label but I'm not able find a way.

To clarify, my problem is the realization of N and S outside the map in the layout with different color:


texts: QGIS DXF measurements in AutoCAD are in Cartesian and not Ellipsoidal

texts:
I'm trying to export some land boundary polygons to AutoCAD using the Export to DXF function.

The issue I am having is that the distance measurements in AutoCAD, of the exported DXF file, represent the same values I get when I select the cartesian option in QGIS measurement and that does not represent the actual land dimensions, I would like to get the ellipsoidal measurements instead, which also match the measurements I get using Google Earth.

I have tried to export the DXF into different CRS, such as EPSG:3857, EPSG:4326 and EPSG:32600 (UTM), but I got the same dimensions.

See images below:


texts: Setting workspace as Python script tool variable/parameter in ArcPy?

texts:
I'm trying to write a script and as I want to automate it as much as possible.

Is it possible to set the workspace as a variable so that when it is run as a standalone script, or as a tool in ArcToolbox, that it will allow everything else to run?

Below, for example, I want the user to select a workspace and then for the script to pick up the shapefiles from this workspace and also create the new geodatabase there.

# Create file Geodatabase
# Import ArcGIS modules
import arcpy
print arcpy.ProductInfo()
# Check out the ArcGIS Spatial Analyst Extension
arcpy.CheckOutExtension("spatial")
# Need to be able to OverWrite Outputs
arcpy.env.overwriteOutput = True
# Set workspace
arcpy.env.workspace = "C:/Prog_Data"
# Set up variables
out_folder_path = "C:/Prog_Data"
out_name = "Prog.gdb"
# Execute CreateFileGDB
arcpy.CreateFileGDB_management(out_folder_path, out_name)
# Set local variables before importing shapefiles to Geodatabase
inFeatures =  ["Points.shp", "Extent.shp", "Rivers.shp", "Population.shp"]
outLocation = "Prog.gdb"
# Execute shapefile to Geodatabase
arcpy.FeatureClassToGeodatabase_conversion(inFeatures, outLocation)


texts: Running ArcGIS Desktop log file in QGIS?

texts: I read that arcgis saves all the processes in a log file, and I was wondering if I can run this log file and all of its processes in a free GIS software like QGIS.


texts: LineString shows only a point at (0.0) of my map

texts:
I want to draw and connect 2 points in my map application. But i get a point at the same location even when i change the coordinates.

addLine(){
var points = [
  [78.65, -32.65], [-98.65, 12.65]
];

let layerLines = new LayerVector({
  source: new SourceVector({
    features: [new Feature({
      geometry: new LineString(points),
      name: 'Line'
    })]
  }),
  style: new Style({
    stroke: new Stroke({
      color: '#ffcc33',
      width: 10
    })
  })
});
this.map.addLayer(layerLines);
console.log(layerLines);
}

The application is developed Angular. you can see the image for more detail of the error.


texts: Python Script for getting elevation difference between two points

texts: I have several stream segments 1000 Km long. I need to find the elevation difference between two consecutive points of distance 1 Km starting from the upstream to downstream. How can I get the elevation difference from DEM? I have stream segments in raster format and also in vector format. It would be better if I got some idea on Python script.


texts: OpenLayers Scaleline always wrong in my virtual star map

texts:
I got a problem for the scale line when I using OpenLayers to show my "star map".

Basically, I give a set of 2D coordinates (X,Y) for each star, and generate the map using MapServer.  Now the map can be correctly shown, but when I add scaleline control, I found the scaleline always starts at 500000000KM which is not correct. The correct scale should be 5400KM, which if I directly use MapServer CGI, it show correctly.

Hope anyone can help me out.

Here is map file:

MAP
    NAME "starmap1"
    PROJECTION
        "init=epsg:4326"
    END
    UNITS meters
    IMAGETYPE PNG
    EXTENT -31230139.51 -8267352.329 36978851.96 9060799.858
    SIZE 800 600`

Here is my JavaScript (using GeoExt):

panel = new GeoExt.MapPanel({
    renderTo: "mappanel",
    map: { maxExtent: new OpenLayers.Bounds(-31230139.510000,-25170977.043091,
        36978851.960000,25964424.572091),
        maxResolution: 85367.949274,
        numZoomLevels:40} ,
        width: 800,
        height: 600,
        floating: true,
        x: 70,
        y: 350
    });

panel.map.zoomToMaxExtent();
panel.map.addControl(new OpenLayers.Control.ScaleLine());


texts: Is it possible to display irregular values in the scale bar in QGIS 2.8?

texts:
Is it possible to display irregular - increasing values in scale bar in print composer?
For example: 5, 20, 100, 250.


texts: Does anybody know a converter or how to convert a geotiff file to a ecw file?


texts: Clipping a raster using an irregular polygon with Python

texts:
I want to clip a raster (GeoTIFF) using a polygon shapefile. Both are in the same reference system. The polygon is not a regular polygon but has many vertices.

I used this Python script below in QGIS 2.18.4

    import processing
    processing.runandload("gdalogr:cliprasterbymasklayer", \
    "/home/giacomo/Desktop/RSM/bordi.tif", \
    "/home/giacomo/Desktop/RSM/mask_vec.shp", \
    "none", \
    False, \
    True, \
    True, \
    "", \
    "", \
    0, \
    0, \
    0, \
    False, \
    "", \
    False, \
    "", \
    "/home/giacomo/Desktop/RSM/bordi_clip.tif")

But I get Error: Wrong parameter value:

The parameters required as input are

    INPUT <ParameterRaster>
    MASK <ParameterVector>
    NO_DATA <ParameterString>
    ALPHA_BAND <ParameterBoolean>
    CROP_TO_CUTLINE <ParameterBoolean>
    KEEP_RESOLUTION <ParameterBoolean>
    RTYPE <ParameterSelection>
    COMPRESS <ParameterSelection>
    JPEGCOMPRESSION <ParameterNumber>
    ZLEVEL <ParameterNumber>
    PREDICTOR <ParameterNumber>
    TILED <ParameterBoolean>
    BIGTIFF <ParameterSelection>
    TFW <ParameterBoolean>
    EXTRA <ParameterString>
    OUTPUT <OutputRaster>


texts: Multiband Raster Reclassification in QGIS

texts: I have a raster which nodata is 0,0,0 in all three bands, how can I reclassify them into 255,255,255?


texts:
Following the documentation for QGIS 3 startup script, I have placed a 'startup.py' script in the following path:

C:\Users\Me\AppData\Roaming\QGIS\QGIS3\profiles\default\python\startup.py

But the script does not get executed. The script only contains the following line:

iface.messageBar().pushMessage("Error", "I'm sorry Dave, I'm afraid I can't do that", level=Qgis.Critical)

When QGIS loads, it should come up with an error message complaing about iface not being recognised. But I receive nothing. I can set up the custom variable and point the path to the Python file but I want to use the startup script method.

Any ideas why?


texts: GeoServer clustering with shared data directory

texts:
I am planning to cluster GeoServer. I can find several results in the web saying that often clustering is set up by configuring multiple GeoServer nodes which all share the same data directory.

At the same time there's some results ( http://osgeo-org.1560.n6.nabble.com/Using-GeoServer-with-shared-data-directory-td4468628.html ) that say that all GeoServer features may not work with shared data directory.

I tried to look at the actual GeoServer documentation but I could not find specifics about how a shared data directory can be used, is it supported, which features work with it and which dont. Best I could find was "When setting up a cluster of GeoServer machines it is common to share a single data directory among all the cluster nodes" at http://docs.geoserver.org/latest/en/user/advanced/logging.html.

It is not a given that a shared data directory can be used. The nodes could be writing all sorts of files in the data dir, blocking and interfering each other. Without clear GeoServer documentation about it I am hesitant to use this approach in production...

Is there some authoritative documentation on

How to use shared data directory with clustered Geoserver
Other aspects of clustering Geoserver (ie. can I simply add a load balancer in front and everything is ok? 3.2 at http://opengeo.org/publications/geoserver-production/ says that, is it correct?)


texts: Re-Projecting lat and long in python Geopandas but geometry unchanged

texts:
Here the question:

the original dataset is like this:

and I want to convert it to normal location like this kind:

I try load the original file in into geopandas and set the projections, to achieve that goal, but unsuccessful.

from geopandas import GeoDataFrame
from shapely.geometry import Point

path_Tw = r'\\Goa\shared\AEP\EMSD\IEAP\EnvData_Prediction\Predictions\Modelling_Projects\Water\Hydrology\1_Data\5_Township\TownshipGrid_Domain\PerGrid(1961-2016)_AllDomain'
path_Tw_loc = os.path.join(path_Tw, "..")

all_tws = pd.read_csv(path_Tw_loc + r'\Township_CentroidPoints.csv', sep=',',index_col=[0])
ft = all_tws.values

point = ft[:,:2]

geometry = [Point(xy) for xy in zip(point[:,0], point[:,1])]
crs = {'init': 'epsg:4326'}
celldf = GeoDataFrame(  crs=crs, geometry=geometry)

point stores the orginal loaded csv information:

I think the geometry will changed to pic2's kind, as I have done this conversion:

geometry = [Point(xy) for xy in zip(point[:,0], point[:,1])]
crs = {'init': 'epsg:4326'}
celldf = GeoDataFrame(  crs=crs, geometry=geometry)

but actually in the final celldf, the geometry still unchanged:

Anyone knows where the problem is?


texts:
I've set up a web service that uses ESRI's REST API spec. This allows me to make requests using the (Silverlight)/WPF SDK.  If I turn on windows authentication though,
I get a 401 error when using GeoProcessor.ExecuteAsynch.

I can successfully call the same url with a WebClient with credentials set:

WebClient c = new WebClient();
c.Credentials = CredentialCache.DefaultCredentials;

If I use WebClient, then I have to deserialize ESRI's json, which is non-trivial, so I'd really like to use Geoprocessor.ExecuteAsynch instead.

Is there some way to send credentials with Geoprocessor?


texts: Error while importing geotiff file into postgresql database in linux platform

texts:
I am trying to import GeoTIFF file into spatially enabled PostgreSQL database.
Using command:

 raster2pgsql -s 4326 -C -F -a path/to/geotiff/file public.table_name
 | psql -d database_name -U postgres.

This commands is working absolutely fine in Windows (Postgres 9.1), but it is giving an error in linux (Redhat 7.2, Postgres 9.1/Postgres9.5).
It gives the following error:

  rt_raster_gdal_warp: could not create GDAL transformation object for output dataset creation
  context: SQL function "ST_transform" statement 1.

I have also tried setting PROJSO and PROJ_LIB environment variable with

  export PROJSO=/path/to/ibproj.so

but it is still not working.


texts: Applying 'less than' filter on list of images in Google Earth Engine

texts:
I am writing a simple GEE code to loop a list of numbers and filter a raster image (DEM) for values less than the numbers in the list.

var aec = function(n){
  var ii = ee.Image(n)    // casting the number to image before filtering but didn't work
  var DEM141 = DEM.lte(ii)

  var DEM141Count = DEM141.reduceRegion({
    geometry: ROI,
    scale: 30,
    reducer: ee.Reducer.sum()
  })
  return DEM141Count
}

var elevs = ee.List.sequence(140,150,1);
var areas = elevs.map(aec)
print(areas)

However, I keep getting this error:
Image.lte, argument 'image2': Invalid type. Expected: Image<unknown bands>. Actual: Float.

The GEE code link is: https://code.earthengine.google.com/6e55360bc82695f930d64ee70a415dc0


texts: Reclassify Raster in ArcGIS Pro 2 by Multiple Bands

texts: I have a raster with three bands (RGB).  I have used the reclassify tool to select only pixels within a certain range of the top band (Red).  Those are assigned 1 and all others are assigned 0.  In order to fine tune my classification, is it possible to use the other bands as well?


texts: Finding nearest neighbors of each coordinate

texts:
I've got a list of coordinates (WGS1984 Decimal Degree) and have started work on a program that performs a DBSCAN on them to find clusters amongst them. To perform the DBSCAN's range query to determine the neighbors within a specified distance, I converted each point into ECEF using formulas from here. I then built a KD Tree to store them.

Is using a KD Tree the best method for this? I am having trouble traversing the tree because points that are clearly within the neighborhood of others are not being found.

Of note, I use the Haversine formula to determine if a point is to be added to a cluster or not, not the square distance between the two ECEF points.


texts:
I wonder if someone knows a way to calculate time-cost isochrones in a bike road network layer, taking the slope into consideration.

The GRASS tool "v.net.iso" allows to create a line-vector layer with categories depending of how far the roads are from a determinated point. Having the data about the speed allowed me to do the analisys with time-cost categories.

But I would like the analisys to show how slope affects the isochrones: It doesn´t take the same time cycling towards the top of a hill than getting down. So taking the slope into consideration would make the analysis more attached to the reality.


texts: Create shapefile containing specific pixels in geotiff

texts: I have a geotiff layer with height data for each pixel (height range = -3.3 to 67.5 m asl). I wish to make a shapefile that contains all pixels within a certain height zone (4.4 to 4.9 m asl).


texts:
I would like to get my print composer map flat for printing and output and look the same as my ArcGIS Pro layout. I have only been using QGIS for a week, trying to replicate my ArcGIS Pro work. Scale in ArcGIS Pro is 1:2,250 and rotation is 120.

Can anyone suggest settings to change?

See pictures of QGIS top and ArcGIS Pro bottom.

I have found that the OSMM topo layers and imagery are in the image coordinate/projection system and most of my other layers are in coordinate system, when a project the features to the same as OSMM the layer ends up miles away, where as in arcgis pro it stays in the same place, I don't want to re-create these features its taken ages.

Any further suggestions


texts: QGIS 2.18.12 REGEX Labelling Dilemma - Carriage Return

texts:
I'm having a drama labelling using RegEx in QGIS.

Ideally, what I want is the lines labelled WITHOUT the 'Spec' in the middle. Which can be a variation of numbers and characters but I need to retain the '1 x' and '24F'

EDIT; 'Spec' can be a variation of characters. e.g 'ABC123' it wont be the word 'Spec'

So the labels will go from this:

2 x Spec - 48F
1 x Spec - 96F
1 x Spec - 12F

To this:

2 x 48F
1 x 96F
1 x 12F

I'm currently using this in the Expression dialog window of the labelling engine;

regexp_replace( "value", (regexp_substr( "value", '[x](.*)?[- ]1')),'')

This works for the multi line values, those with a return. But it eliminates anything on a single line.

I need the single line values as well as the multi line strings.

Any pointers on what I could be missing?


texts:
For complex labels, it would be useful to be able to use HTML formatting when different fields should be shown, on multiple lines and with different font sizes.
Is this possible with QGIS? Or is there a workaround?


texts:
I am working with QGIS and I have a DTM of a river. I need to extract multiple cross sections (one every 10m), and they must be perpendicular to the river center line.
Can you tell me how to do that, with the passages?


texts: using pycuda to accelerate the processes work in Qgis

texts:
I wrote a module in QGIS which I work with the big raster .tiff layers.

I want it to use the video card for processing (not the CPU).

Are there are any ideas how "processing.runalg" can be run in parallel on the GPU?

 def change(self):
        var=self.dlg.lineEdit_3.text()
        var2 = self.dlg.lineEdit_5.text()
        processing.runalg('saga:changegridvalues',var, 1, "-100,0,0",var2)
        tif4= self.dlg.lineEdit_5.text()
        layer = QgsRasterLayer(tif4, "rastr_change1")

For example using PyCUDA.


texts:
After hours and hours of reading I was finally able to get multi-line labels established on my current project. I'm using QGIS 3.16.2 and I'm trying to get multi-line labels positioned at different positions around an item. As you can see in the image below these are a few of the labels I want to use and I'd like for them to be at different positions around the well symbol. For example the long string of numbers I want posted above, for the operator name (second label) I want it to the right of the symbol, the well name (third label) below the symbol, and lastly the well number (fourth label) to be to the left of the symbol.
I've read and read and I can't seem to find anywhere where this has been discussed or solved.


texts: Render vectors from OSM data in POSTGIS database

texts: I'm a newbie when it comes to mapping application. I had successfully loaded my OpenstreetMap data in a Postgresql/PostGIS database. Now my problem is how to render the vector data (the POI's, ways, etc) from the database to the the map. I don't have the slightest idea on how to do it, though I think I knew how they are rendered (OpenLayers use KML/XML files to render vector data, I tried it in their tutorial). Can someone explain to me how it works; do I need to get the data from the database as XML file? I don't know if I'm missing here something. Please help me. Thanks a lot!


texts:
As ArcGIS 10 doesn't really support the full COM, I am looking at creating an extension.

I am looking for samples to look through, to familiarise myself with them, but cannot find any.

Any chance someone can show some links ot some samples of creating extensions

Using VS2010 and C#


texts:
I am processing sentinel 2 data and while most bands look fine I have noticed an odd seam in band 1 of several images as can be seen from the attached image.

The image is cut from the ESA scene S2A_MSIL2A_20190921T101031_N0213_R022_T33VWD_20190921T130515 and displayed using only the value range 0-300 of band 1.

What type of artefact is this seam? Is it an error with my processing or a property of the instrument or what is it?


texts: Tooltips are not captured by leaflet-image

texts:
I am trying to capture my map with leaflet-image, but the tooltips are not captured. I enabled preferCanvas in the mapOptions and everything except from the tooltips is captured. Each one of my markers has tooltips. The tooltips seem to be drawn in HTML. I have even tried to do something like this:

var canvas = document.createElement("canvas");
var ctx = canvas.getContext("2d");
ctx.font = "30px Arial";
ctx.fillText(feature.properties.test_number,10,50);
return canvas;

The numbers are now in canvas. But the wrapper still is in HTML - I believe by default from Leaflet. I have no problem modifying my tooltips (they're just numbers) to be in canvas or modify the way I am doing things in general to make the capture possible.

Why are tooltips not being captured in screenshot?

My goal is the  ability for leaflet-image to be able to capture the tooltips that are bound to the markers.

Map Doesn't render (due to Bing Maps issue -- different Issue)


texts: Improving performance when importing file geodatabase into PostgreSQL?

texts:
I have a file geodatabase created with ArcMap 10.1 that contains one feature class with 650 million points. The feature class contains the shape field and an identifier. The .gdb is approximately 26GB.

I'm using ogr2ogr to import the .gdb into a PostgreSQL database. I started the process and the features are being inserted into the database. I verified this with a SELECT COUNT(*) FROM <TABLE> and the number of rows is increasing.

Based on my rudimentary approach to timing - watch the PC clock and use the above SELECT - I estimate the import is progressing at 7 million features per hour. Some quick mental math and I'm looking at 90+ hours to import the .gdb - if the rate remains the same and doesn't degrade over time. This is roughly 2000 features per second.

The database is an out of the box PostgreSQL 9.4.5 / PostGIS 2.2 installation on a Windows 10 PC with a SSD (400 GB free / 512GB total) and 16GB RAM. No PostgreSQL configuration changes have been made yet as I'm not sure what, if anything, I should do.

Is there anything I can do to increase performance of the .gdb import? Perhaps a recommendation that maximizes writes and  then I can rollback the changes after the imports? I could let this process run until completion but there will be other imports I need to perform and I'd rather not wait four days for each one. I have full control of this PC so I'll cautiously say I'm open to any suggestions.

Update:
I'm using the following command and options. I don't know how to check if this is the open or ESRI library. The version of GDAL is 1.12 (95% sure) and was installed with QGIS 2.12 Lyon (100% sure of this).

ogr2ogr -f "PostgreSQL" PG:"host=localhost port=5432 dbname=spatial_playground user=postgres password=********" c:\gis_data\pointdata\points.gdb

The database and .gdb are on the same hard drive so contention is probably impacting this. I considered throwing the .gdb on a USB 3.0 hard drive. The OGR command is running on the same machine. Everything is being done on the same machine; No network involved yet.


texts: Overwrite ArcGIS geocode service after editing the attribute table on arcmap

texts: i have an ArcGIS geocode service and it's running successfully. Now I need to edit the attribute tables to add new rows and overwrite the service. how can I do this?


texts:
I have ArcGIS 10.5 and some add-ins which I want to be automatically updated after making some corrections. Recently I have found a function that provides ArcGIS to do it. This is all about Add-in Manager - setting custom folder with add-ins.
Once I set a folder and my tools appeared but now I want to change or delete this folder. I removed this folder from box in Add-In Manager and after restarting ArcMap nothing changed. The path to a folder still exists despite I removed it manually.
Almost the same thing about other users - they are not able to save this path and as a result they are not able to setup tools automatically.
Looks like the thing is about overwriting some of ArcGIS parameters, it blocks them. On 10.2 everything works pretty well.

Is there any solution?


texts:
How to create markers/icons at the startpoint and endpoint automatically in OL/JS? I'm doing this right now - I know startpoint and endpoint from PHP analysis:

var startpoint = new ol.Feature({
  geometry: new ol.geom.Point(
    ol.proj.fromLonLat([xx.xxx,xx.xxx])
  )
});
var vectorSource = new ol.source.Vector({
  features: [startpoint]
});
var start = new ol.layer.Vector({
  source: vectorSource,
  renderBuffer: 0
});
start.setStyle(new ol.style.Style({
            image: new ol.style.Icon({
            color: 'green',
            src: '/gpx/start.png',
            anchor: [0.6, 1]        })
        }));
map.addLayer(start);

So I'm trying something like this:

source.on("drawstart", function(e1) {
var startpoint=map.getEventCoordinate(e1);`

...

start.setStyle(new ol.style.Style({
            image: new ol.style.Icon({
            color: 'green',
            src: '/gpx/start.png',
            anchor: [0.6, 1]        })
        }));
map.addLayer(start);
});

but without success.