texts: How to handle display of nested polygons (census, county/state, etc) for a given map extent

texts:
My question is how best to decide which layer (from a set of nested, aggregated census polygons) to display for any given map extent in an online map.

Australia's census boundaries are built from multiple layers, aggregated from highly detailed collection districts through a number of intermediate layers to eventually comprise the whole country. This same principle applies to city/county/state/country boundaries.

Switching census layers based on scale dependency is unsatisfactory, since there must be a compromise between the populated cities and the vastly empty interior. Eg, over Sydney the Collection Districts should show at around 1:10k whereas in rural areas they should show above 1:500k

1) Are there any examples available in online maps where this problem has been solved elegantly? All of the examples I've seen use scale dependencies to switch layers on/off

2) I'm contemplating counting the features every time the map extent changes, and using this to decide which is the best layer to display. This "breaks" the traditional GIS paradigm whereby layers can switch on/off based on scale, since in my scenario they may also switch on/off based on different map extents at the same scale.

Are there any examples of this in action? Does anyone have any thoughts or advice on whether this will be confusing to the user?

Eg, a user might be viewing Collection Districts, then pan the map slightly. This could cause the Collection Districts to switch off, and the next layer to switch on. Has this been done anywhere before?

Thanks,
Steve

I've put a sample map at http://203.14.35.213/atlas-2.0/public/nsw/home/map/australia-population.html showing the direction I'm hoping to go. Note that this is a QA server so it may be slow, and not up 100% of the time.

Pan and zoom around the map down to street level, and note that the polygon boundaries change (if you open the Legend panel, you can see which of the 5 layers is currently displayed).

The boundary layer is chosen each time the map draws, based on the count of features rather than the scale (eg, a different boundary will draw over Sydney, compared to a rural area at the same scale). This can mean that panning the map results in the layer changing.

The idea is for this process to be transparent to the (novice) user, who should just see the "best" information wherever they browse.

My original question was whether the experts feel this "breaks" the traditional model of scale-dependent rendering. Does anyone have any feedback on whether this approach has merit or is flawed?


texts: Workflow for generation of contours with LiDAR data and OCAD

texts:
I am updating our maps for an orienteering club in Illinois. I was very excited to find out that we have access to .las lidar data in 2000x2000ft files. I checked out OCAD11 and it can handle .las files so I thought I was golden.  Except the .las files are in feet and OCAD11 does not handle feet and the conversion to utm is not right. I have contacted OCAD and they are unable to solve this issue. So I am looking for a possible work around or alternate workflow for develop georeference contours in utm.

I have looked at Fusion and also the article "Using Lidar point clouds for Orienteering base map generation" by Terje Mathisen. I do not see a clean path for dealing with multiple .las files covering a park and also dealing with the Illinois state plane to utm.

My current plan is to create a dem in OCAD11 export to a dem file. Import the file to QGIS, do the conversation to utm, export to xyz file, import back to  OCAD, create contours.

Can someone suggest a cleaner simple workflow? I would like to avoid buying OCAD11. We currently use OCAD9. I would like to be able to make the contours externally and import them into OCAD9.


texts: Display overlapping polyline segments as multiple lines with offsets?

texts:
Selecting which polyline is displayed on top in case of segment overlaps seems to be quite a FAQ, however, I would like to find a way to display overlapping polyline segments side-by-side, i.e. as multiple lines separated by offsets, similar to how Google Transit displays branching transit services.

Is there a Leaflet library that does that?


texts: Adds ArcGIS controls in the Visual Studio Toolbox window for dragging and dropping into a design view surface

texts:
I want to drag and drop MapControl into a design view, but I have an error: System.Runtime.InteropeServices.COMException(0x80040154):Classe non enregistrée (Exception from HRESULT:80040154(REGDB_E_CLASSNOTREG)).
 I find that your map control is not registered or Class not registered,
but I don't know how to do it
in VS 2008 SDK 9.3, please help.


texts:
I have a local ENU coordinate system on the surface of the Earth (similar to this question: GPS coordinates to XYZ on relative grid). I know the lat,lon,alt of the local ENU origin, and thus I can also get the origin in terms of ECEF.

Now, I have a point of interest that I'd like to know its ECEF XYZ coordinates. I don't know the lat/lon/alt of the POI, but I do know its position relative to my local ENU grid/origin. How can I calculate the ECEF XYZ coordinates of my POI? An alternative approach might be to first get the lat/lon/alt of the POI.

I tried using this formula (https://gssc.esa.int/navipedia/index.php/Transformations_between_ECEF_and_ENU_coordinates) to convert ENU to ECEF, but it's not really working. It's obvious that the formula has no consideration for the fact that the ENU origin has an altitude. I tried applying some earth radius distance to the "U", and that can get me close to the correct answer, but it looks like this also changes the latitude, so I assume I'm using an incorrect approach.


texts: Sort data in an ascending order based on the geometry of an sf object in R (Sort by Lat/Long)

texts:
I am trying to sort my data based on the latitude and longitude. I have an sf object in R, which essentially stores the latitude and longitude as a separate column called geometry.

See below:

> head(barfly.pr.sp)
Simple feature collection with 6 features and 14 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 76.99832 ymin: 10.04038 xmax: 77.0799 ymax: 10.1619
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

    effort_distance_km effort_area_ha reviewed                  geometry
0                 0.00              0    FALSE  POINT (77.0351 10.14666)
0.1              24.14              0    FALSE   POINT (77.066 10.08324)
0.2               8.00              0    FALSE POINT (77.05931 10.09309)
0.3               3.00              0    FALSE  POINT (77.0799 10.04038)
0.4               1.20              0    FALSE  POINT (76.99832 10.1619)
0.5               2.50              0    FALSE POINT (77.06757 10.15112)

I am trying to sort the same based on geometry and I get this error:

> barfly.pr.sp <-  barfly.pr.sp[order(barfly.pr.sp$geometry),]
Error in Ops.sfc(xi, xj) : operation > not supported

I also tried dplyr::group_by() and I got the following error:

barfly.pr.sp <-  barfly.pr.sp %>% group_by(geometry)
Error in grouped_df_impl(data, unname(vars), drop) :
  Column `geometry` can't be used as a grouping variable because it's a sfc_POINT/sfc

Any suggestions?


texts: I have a georeferenced raster with water depth values and I would like to add more pixels in order to extend my raster as far as the coastline (like in the picture, dark blue has no values). How can I do that?


texts: Convert Sql Server Spatial data to GeoJSON & plot on map

texts:
I have a sql server table with Geo spatial data. An example of row would be (with Column names):

type: streetline
code: 231001
geog: 0xE6100000011 ........
Centroid 0xE61000000C.......
geom: 0xE6100000011 ........

I am using C# with MVC3 to draw the above shapes on leaflet map. I am fetching the above sql data into a datatable. Using the sql function "ToString()" on column "geom" I get the following:

"LINESTRING (-1.131510412 52.65531, -1.13286 52.65559)",
"POLYGON ((-1.1116360 52.6409953, -1.1116683 52.6413, -1.11146723 52.641317, -1.11133263 52.6413572, -1.1113059))",

The question is how do I convert the above to GeoJSON to be able to plot on the leaflet map. Here is an example of the expected output (GeoJSON) I want:

var geojsonFeature = {
    "type": "Feature",
    "properties": {
        "name": "Coors Field",
        "amenity": "Baseball Stadium",
        "popupContent": "This is where the Rockies play!"
    },
    "geometry": {
        "type": "Point",
        "coordinates": [-104.99404, 39.75621]
    }
};


texts:
I'd like to convert the big red area to a smaller one (blue) by shifting their borders parallelly by a certain distance to inner direction. This is not equivalent to scale as the result shape will look different.
Do you know an existing out-of-the-box function for contract? I tried to find one e.g. on turf but I did not manage to.


texts:
I am trying to read a GDAL VRT file.

I have included

    <!-- VRT Format (reader) -->
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-imageio-ext-gdal</artifactId>
        <version>${geotools.version}</version>
    </dependency>

and I have these files:

VRTFormat
VRTFormatFactory
VRTReader

https://github.com/geotools/geotools/tree/669897d2778dab8c31e7ad97bf16ae02e885bfa1/modules/plugin/imageio-ext-gdal/src/main/java/org/geotools/coverageio/gdal/vrt

but

 AbstractGridFormat format = GridFormatFinder.findFormat("...");
 GridCoverage2DReader reader = format.getReader("...");

causes Exception

Trying to get a reader from an unknown format.
java.lang.UnsupportedOperationException: Trying to get a reader from an unknown format.
    at org.geotools.coverage.grid.io.UnknownFormat.getReader(UnknownFormat.java:62)

how to enable GeoTools to find the respective VRT reader?


texts: ESRI Flex API - Dealing With Graphics Layers in Polar Regions

texts:
I have developed a Flex app that incorporates a number of geo-processing functions.  The Base Layer of this app is the standard ESRI World Imagery in Web Mercator coordinate space.

One of these geo-processing functions is a glorified buffer service that can (and often) gets run in or near the polar regions.  The graphics in these locations run-off the map (see attached image)

So, I have 2 questions about dealing with this situation:

Is there some way to confine the graphics layer to the map surface,
just like ArcMap does?
Is there an algorithm I could use (client side) to draw the graphics
like they would fall on an actual globe?


texts: EPSG codes for projections to use when computing the area of global polygon sets in python env

texts:
The general answer to this question is "Use any equal area projection in units of meters that covers the world". Python code needed to make use of this is short and sweet:

import geopandas as gpd

gdf = gpd.read_file("my.geojson")
equal_area_crs = "EPSG:????"
areas = gdf.to_crs(equal_area_crs).area

I found myself very frustrated however that googling this phrase (and many variants of it) produced zero useful results on the first page. There are MANY introductions to the basics of projections but no quick answers I could use to quickly fix a bug that used an improper area distorting projection.

Searching on epsg.io is somewhat inflexible, seeming only to search in names. A search for "equal area" or "lambert azimuthal equal area" and a variety of other equal area searches produce results with specific regional applicability and mixed units.

I'm posting here because it took me a full hour to find a satisfactory answer, and this may be useful to others.


texts: Using Arcade to label, using coded values in ArcGIS Pro

texts:
I am trying to label features that have coded domains in ArcGIS Pro. In ArcMap the procedure is to use the coded domains in an expression as shown below, but you have to enable the codes using the checkbox; as shown in the snippet below. In ArcGIS Pro this is Arcade driven and we are having issues coming up with a correct expression to achieve this.

How would I do this using Arcade in ArcGIS Pro?


texts:
I gathered field points with a TOPCON hand held receiver.  I saved the points as a shapefile and brought them into ArcMap.  I used the editor tool to connect the points and created a polygon.

Is there a tool that I can use within ArcMap to find the perimeter of the polygon?


texts: Issue Trying to create Zonal Statistics using Gdal and Python.

texts:
I am trying to create zonal statistics using Python and Gdal. I have a polygon shapefile and a raster file, and in order to do so, I am using a piece of code I found in StackExchange.

The raster and shapefile used can be found here

Here is the code I am using:

#!/usr/bin/python
#coding: utf-8

## Code from Stack Exchange: http://gis.stackexchange.com/questions/21888/how-to-overlay-shapefile-and-raster
## User: ustroetz

# Calculates statistics (mean) on values of a raster within the zones of an polygon shapefile

import gdal, ogr, osr, numpy

def zonal_stats(input_value_raster, input_zone_polygon):

 # Open data
    raster = gdal.Open(input_value_raster)
    driver = ogr.GetDriverByName('ESRI Shapefile')
    shp = driver.Open(input_zone_polygon)
    lyr = shp.GetLayer()

    # get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # reproject geometry to same projection as raster
    sourceSR = lyr.GetSpatialRef()
    targetSR = osr.SpatialReference()
    targetSR.ImportFromWkt(raster.GetProjectionRef())
    coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()
    geom.Transform(coordTrans)

    # Get extent of geometry
    ring = geom.GetGeometryRef(0)
    numpoints = ring.GetPointCount()
    pointsX = []; pointsY = []
    for p in range(numpoints):
            lon, lat, z = ring.GetPoint(p)
            pointsX.append(lon)
            pointsY.append(lat)
    xmin = min(pointsX)
    xmax = max(pointsX)
    ymin = min(pointsY)
    ymax = max(pointsY)

    # Specify offset and rows and columns to read
    xoff = int((xmin - xOrigin)/pixelWidth)
    yoff = int((yOrigin - ymax)/pixelWidth)
    xcount = int((xmax - xmin)/pixelWidth)+1
    ycount = int((ymax - ymin)/pixelWidth)+1

    # create memory target raster
    target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, gdal.GDT_Byte)
    target_ds.SetGeoTransform((
        xmin, pixelWidth, 0,
        ymax, 0, pixelHeight,
    ))

    # create for target raster the same projection as for the value raster
    raster_srs = osr.SpatialReference()
    raster_srs.ImportFromWkt(raster.GetProjectionRef())
    target_ds.SetProjection(raster_srs.ExportToWkt())

    # rasterize zone polygon to raster
    gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])

    # read raster as arrays
    banddataraster = raster.GetRasterBand(1)
    dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float)

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float)

    # mask zone of raster
    zoneraster = numpy.ma.masked_array(dataraster,  numpy.logical_not(datamask))

    # calculate mean of zonal raster
    return numpy.mean(zoneraster)

However, when I implement it, the output (numpy.mean(zoneraster)) is just a number and not an array or a raster. What am I missing? Is there something wrong with the files imput? The code?


texts: How to install pgrouting from source on CentOS 6.3?

texts:
I'm running into a problem once I try

[cmake_dir]/bin/cmake .

I get these errors:

Boost headers were found here: /usr/include
Output directory for libraries is set to sh: /usr/bin/postgres/pg_config: Not a directory
Installation directory for libraries is set to sh: /usr/bin/postgres/pg_config: Not a directory and for SQL files is set to /usr/share/pgrouting
-- Configuring done
-- Generating done

I'm running CentOS 6.3 and all I know is that my postgresql data folder is at /var/lib/pgsql/data


texts:
Possible Duplicate:
  Cloudmade Web Maps API  Marker Clustering problem

I'm at Geocoding in Cloudmade Web Maps API. I've the code below. Whenever I try to access the page, content of elements must consist of well-formed character... trace is generated, and the page does not render.

My environment is Apache Tomcat 7, and JSF.

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

  <title>CloudMade Geocoding</title>

I've, of course, replaced the API_KEY with my actual API key, which works with the rest of the examples. Is anyone able to spot what I may be doing wrong.

Thanks in advance.


texts: Plotting lon/lat point (without manual conversion) on XY map using GeoPandas?

texts:
I am data scientist not experienced in GIS and need to do some simple map plotting of the Netherlands. Using Geopandas everything works out nice and dandy for my purposes however I am stuck on 1 aspect:

The map I drew has the following CRS:

map_df__.crs
{'init': 'epsg:28992'} (this is what NL looks like to de dutch native)

I want to plot a couple of points of which I get the lon/lat coordinates from Google Maps.

Muddling my way to understand projections on the fly I got to this code:

map_df.to_crs(epsg=4289).plot(edgecolor='black')

Which gives this map: (Red dot: Amsterdam is an example)

To a Dutch person this map looks a bit squashed. I am guessing because because of the projection.

How do I plot a lon/lat point (without manual conversion) on an XY map i.e. how can make the squashed map look like the 'stretched' version?

The suggested espg:
map_df.to_crs(epsg=3857).plot(edgecolor='black'), gives:


texts:
I'm new to QGIS.

I have imported a shapefile for a road network. I have then imported a table containing condition data about the network. I have used the LRS plugin and the data displays perfectly and the "Attribute table" is full of data.  I have saved the project, but then when I close QGIS, then reopen it and open the project, whilst the condition layers are listed in the layers panel, none of its data is displayed and its "attribute table" is empty.

How do I save the project so that when I open it, I can continue where I left of rather than having to recreate everything?


texts: Combine multiple polygons into a single one (with holes)

texts:
I have the following GeoDataFrame with multiple polygons which represents a single area (raster_val = 1.0) with holes in the middle (raster_val = NaN):

                                       geometry           raster_val
0   POLYGON ((8.77500 47.70000, 8.77500 47.67500, ...         NaN
1   POLYGON ((8.90000 47.40000, 8.90000 47.30000, ...         NaN
...
9   POLYGON ((6.12500 46.27500, 6.12500 46.25000, ...         1.0
10  POLYGON ((8.57500 47.80000, 8.57500 47.77500, ...         1.0
11  POLYGON ((8.90000 45.85000, 8.90000 45.82500, ...         1.0

I would like to combine polygons corresponding to raster_val = 1.0 into a single Polygon/Multipolygon but preserving the holes indicated by the polygons with raster_val = NaN.

I've already tried geopandas dissolve() method:

polygons = polygons.dissolve(by='raster_val')

This combines raster_val = 1.0 polygons but it ignores all raster_val = NaN polygons so the combined multipolygon doesn't have any holes. Do you have any suggestions on how to include those NaN holes?


texts:
I have a system where my users will be uploading datasets that include a number of attributes proprietary to my system but may also include street addresses. At the time of upload, I would like to geocode these addresses and store the coordinates in my database. This is needed so that at a later time, other users can visualize the address data on a map and draw a "geofence" around a certain area to perform actions on only a subset of this data.
The trouble I'm having is that the major geocoding providers (Google, Mapbox, HERE etc) seem to prohibit this usage because it requires storing the results of geocodes for a long time. And I can't very well re-geocode all the addresses on the fly when a user wants to see some map data - that could be thousands of records at once!
How does one handle a case like this? Do I need Mapbox's "permanent geocodes" or whatever the other vendors call them? Or is there some other way?
Note - most, but definitely not all, of my customers who need this will be located in the USA.


texts: shp to dgn microstation export file error

texts:
I want to export my shapefile into microstation dgn using qgis. But when I do click right on layer and "save as" this message error shows: Export to vector file failed. OGR ORIGIN is not a valid 2d or 3d tuple. Separate tuple values with comma.

Am I doing something wrong?


texts:
I want to plot a monthly SO2 concentration for a city day by day in Google Earth Engine. I can visualize on the map by layer but I can't plot a graph.

This is the code I have written:

var antalya=turkiye.filterMetadata('NAME_1','equals','Antalya');
Map.addLayer(antalya);

var before_covid_start=('2019-04-01');
var before_covid_end=('2019-04-30');
var after_covid_start=('2020-04-01');
var after_covid_end=('2020-04-30');

var SO2_concentration=ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_SO2')
        .filterBounds(antalya)
        .select('SO2_column_number_density');

var before_SO2 = SO2_concentration.filterDate(before_covid_start,before_covid_end);
var after_SO2 = SO2_concentration.filterDate(after_covid_start, after_covid_end);

var beforeSO2=before_SO2.median().clip(antalya);
var afterSO2=after_SO2.median().clip(antalya);

Map.addLayer(beforeSO2,{
  max:0.0005,
  min:0,
  palette: ['black', 'purple', 'blue', 'green','yellow','orange','red']}, 'SO2 Total - 2019', true,0.6);

Map.addLayer(afterSO2,{
  max:0.0005,
  min:0,
  palette: ['black', 'purple', 'blue', 'green','yellow','orange','red']}, 'SO2 Total - 2020', true,0.6);


texts: What is the difference between spatial regimes and spatial heterogeneity?

texts:
Let me cite excerpts from the GeoDA-Gloassary (https://geodacenter.asu.edu/node/390)

Spatial regimes:

  Spatial regimes are a form of spatial heterogeneity, which implies structural differences across space. When a variable is characterized by distinct distributions (e.g., with a different mean or variance) for different geographic subregions, these subregions might point to the existence of spatial regimes.

Spatial heterogeneity:

  Spatial heterogeneity exists when structural changes related to location exist in a dataset. In such cases, spatial regimes might be present, which are characterized by differing parameter values or functional forms.

Unfortunately, the difference is still not clear to me. Both explanations seem to indicate that spatial regimes as well as spatial heterogeneity refer to a process' variation of mean/variance over space.


texts: Using Python interpreter other than that of QGIS interpreter?

texts:
This is a straightforward question but I've been looking for a solution for such a long time. The document "IDE Settings for Developing QGIS" says: 'on Windows you need to make sure you that you have the same environment settings and use the same libraries and interpreter as QGIS.'

How can I use another Python interpreter (e.g the one in some Python distributions such as Winpython) with the QGIS environment settings (I almost know how to do that by using batch files to set the correct environment settings) and then how can I use the QGIS essential libraries?

I'm asking this for the Windows operating system. If anyone has any idea on how to do it for QGIS 3.x that is much more preferred.


texts: QGIS CRS correct, but layers appear warped

texts:
I am trying to overlay a raster Landsat 8 scene and vector scene boundary. Both are using WGS 84 + UTM48N. However, the vector (orange line) appears offset when overlaid on the raster. I have checked and confirmed that both layers have the same CRS (EPSG:32648), which is also being used by the QGIS map project.

A similar issue has been brought up before, e.g. here, which noted a bug, but nearly two years later I assume this is now solved.

Anyone else experiencing this and can confirm if it is a bug or a user error?


texts: I have a problem. My objective is to do the cache of a WMS service, imagining that I have a layer with a digit ('1' or '2'), and have a style that says for example that '1' is the color blue and '2' is green. My problem is that if I change for example the feature from 1 to 2, the cached files are still in blue and don't update to green. How can I force the erase of the cached files if the features of the service changes?


texts:
I have a drone image is which is georeferenced over a Google Sat. image.   These align well, and both use WGS84-Pseudo Mercator (EPSG:3857).

I have a set of data points in a shapefile which are provided as WGS-84 (EPSG:4326); and QGIS is set for on-the-fly re-projection.   Overlaying the data points over the images, I notice that all the data points are offset from the underlying physical objects seen on the images -- for the most part the shift is consistent across the entire data set of 20 or so points (allowing for GPS errors at the time of data collection, I suppose)

My question is this:

How does one physically shift the data set, that is by shifting the CRS to a custom CRS, by that small offset?

I would assume that I create a custom CRS, based on the existing WGS84, and then add some terms to the command line which accomplish the shift in the X and Y directions, by the 10-20 feet.  However, I can't find an instruction tool/set online that explains this.  Help...

Alternatively, I would have hoped for a way to merely enter a vector to define that shift (that is, to create a vector from a data point to the what should be the underlying physical object in the drone photo)... but I don't see that capacity in the software.

Yet another alternative is that there is some odd/weird CRS in the data collection which resolves the shift, and/or a weirdness in the WGS-84 pseudo-Mercator projection around where my data lives on the ground. So, I wonder if there is a way to discover which other formally named projections in QGIS live within 10 feet of my data points (and therefore by using that existing alternate, I get my data points shifted).

Of course, I can't disallow that the underlying Satellite image is, as a base map, shifted.   In my use, in this situation, I could shift the GE Image, using a similar CRS adjustment.  So... same questions.

I am using QGIS 2.6.1.  I realize I could drag each individual data point, manually.  Part learning curve, with this question, in advance of a task that will have far more points.


texts: Leaflet/JS: Adding polylines after pan/zoom

texts:
I am using a Leaflet map as an image viewer, and I'm trying to add a 'toggle crosshairs' function that will place a crosshair at the center of the image. I've got it working to a point where the crosshair appears correctly, and is 'coupled' to the image when doing pans/zooms, which is what I would expect to see.

The issue I'm having is that if I do the pan/zoom FIRST, and then toggle the crosshairs, the crosshairs appear in the center of the 'window', not in the center of the image.

The code I'm using for this is here:

    toggleCrosshair(toggleVal) {
    let map = this.map;

    if (toggleVal == true) {
        var north = map.getBounds().getNorth(); // these are in lat, lng
        var south = map.getBounds().getSouth();
        var east = map.getBounds().getEast();
        var west = map.getBounds().getWest();

        var center = map.getCenter();

        var latlngHorizontalLeft = new L.latLng(center.lat, west);
        var latlngHorizontalRight = new L.latLng(center.lat, east);

        var horizLineList = [latlngHorizontalLeft, latlngHorizontalRight];
        let horizLine = new L.Polyline(horizLineList, {
            color: 'red',
            weight: 2,
            opacity: 0.5,
            smoothFactor: 1
        });

        var latlngVerticalTop = new L.latLng(north, center.lng);
        var latlngVerticalBottom = new L.latLng(south, center.lng);

        var vertLineList = [latlngVerticalTop, latlngVerticalBottom];
        let vertLine = new L.Polyline(vertLineList, {
            color: '#A5CE3A',
            weight: 2,
            opacity: 0.5,
            smoothFactor: 1
        });

        horizLine.addTo(map);
        vertLine.addTo(map);

    } else {
        var i;
        for (i in map._layers) {
            if (map._layers[i].options.color == '#A5CE3A') {
                try {
                    map.removeLayer(map._layers[i]);
                } catch (e) {
                    console.log("problem with " + e + map._layers[i]);
                }
            } else if (map._layers[i].options.color == 'red') {
                try {
                    map.removeLayer(map._layers[i]);
                } catch (e) {
                    console.log("problem with " + e + map._layers[i]);
                }
            }
        }
    }
}

I've asked this question StackOverflow, but haven't gotten a response, so I figured I would try here.


texts: QGIS print composer internal grid labels

texts:
Is there a way to have the grid labelled inside of the map frame (I don't mean the option to draw coordinates inside the frame)? What I want to achieve is a "ladder" of labels like in this example:

http://www.freytagberndt.com/wp-content/uploads/2015/12/9007868048732_Blattschnitt03.jpg

Not only near the borders of the frame but also repeatedly inside (because the maps are usually quite large). That's why the "Inside frame" option is not sufficient.


texts: How to create animated timeline in ArcGIS for Desktop?

texts:
I can't find any options in ArcGIS timeslider to create a timeline like presented in this video. The video itself is a guide for creating time animation, but no information is given how to create such timeline like in the end.

I have tried to search some information about that and can't find anything usefull.

Does ArcGIS support such functionality?

UPDATE 1:
I need an autorefreshed timeline (which is presented not as usual date text).
It should be animated line where each date is a vertical line on a timeline.
Also I need to label each date on it as a text. See picture below (frame by frame animation).

Added a post at ArcGIS Idea, to add such functionality. Upvotes are welcome.


texts:
I downloaded a file .osm. I want to have the same file but with extension .xml
To do that, I downloaded osmconverter but the problem is, as states in wikipedia I can convert the .osm file into specific formates and .xml is not one of them.

Please let me know how to do that because I need to do further processing on the file and it must be .xml


texts: Multiple Reclassifications by a variable

texts: I am using Model Builder to reclassify a DEM raster (the input raster is a parameter and may vary in range of elevations). (The model then goes on to extract the highest point in each patch above a given elevation.) I can classify above and below a fixed height, but what I would like to do is make multiple passes, incrementing the classification cutoff by, say, 10 each pass (or pass the increment as a parameter). How can I change the Old Values line to accept the changing cutoff? I see there is a Load and Save function, but I am unclear how I would write a table within Model Builder. Ideally the model would have the start elevation, the interval, and the final elevation as input parameters. I can do this in Python, but would like to have it in Model Builder format.


texts:
I'm trying to get features from WFS service, this service using stored queries. I've managed to retrieve data both with axios and XMLHttpRequest. But the problem in the reading features and showing them on my layer.

My sample code:

let vector;

let xmlhttp = new XMLHttpRequest();
xmlhttp.onload = function() {
    var format = new GML3();
    var xmlDoc = xmlhttp.responseXML;
    vector = new Vector({
        source: new VectorSource({
            format: format
        })
    });

    for (var i = 1; i < xmlDoc.children[0].children.length; i++) {

        console.log('xmlDoc.children[0].children[i] >>> ', xmlDoc.children[0].children[i]);
        let features = format.readFeatures(xmlDoc.children[0].children[i]);

        console.log('features >>> ', features);
        // features.getGeometry().transform('EPSG:4326', 'EPSG:3857');
        vector.getSource().addFeature(features);
    }

    map.addLayer(vector);
    // map.getView().fit(vector.getSource().getExtent(), map.getSize())
};

xmlhttp.open("GET",url,true);

// trigger the GET request
xmlhttp.send();

So when i run this code, console tells that "Uncaught TypeError: feature.getId is not a function". I'm using OpenLayers 4. Could anyone explain what happening ? As i understood there is a class Feature , which has set and get id's methods. Moreover is it possible to convert my GML response into the geoJson , maybe it will be easier to work with ?


texts: Where can I see a list of all renamed functions?

texts:
I'm clearing a large project of all legacy PostGIS code, and I'm looking for a comprehensive list of everything that could be a deprecated function call. I'm not going to change it all in one swoop -- I'm aware that POINT() is probably in some WKT and should not be replaced with ST_POINT(). But I'd still like to know them.

All the pages I've found so far just say "Stop using the deprecated functions!" and lists the new ones but I don't know for sure what every single one of the old ones is. It's not as simple as "All new ones without the prefix" because, for instance, I doubt there could have been a PostGIS COUNT() that was improved to ST_COUNT().


texts: How to tell zonal statistics plugin to ignore "No Data"?

texts:
I have a shapefile and a raster. They superimpose to a greater degree except at edges where raster pixels do not extend up to layer boundary (see image if it was not clear).  The green line is from the shapefile. The red color is raster. The space between them is a NO DATA.

It must probably because of Larger pixel sizes of rasters that caused problem when I clipped the raster.  The problem arises when I use zonal statistics. This tool calculated mean considering the NO DATA field values as well.

  Is there a way to tell this plugin (or may be another better plugin) to ignore No Data from the raster while calculating mean?


texts: Created GeoTIFF only has NaN values in PostgresSQL

texts:
I created a GeoTIFF file from a numpy array, loaded it into a PostgreSQL DB and now I want to retrieve this raster data and convert it back into a numpy array. I know this seems nothing really practically, but its just exemplary for another bigger task.
But when I create the numpy file, it only loads NaN Values. I checked creating the numpy array with another GeoTIFF file (as they describe here: https://gist.github.com/gzagatti/a649674d85b480b98b7d75b31d4c9555)
So the mistake must be happened, when I created the GeoTIFF.
# Create dataset file
drv = gdal.GetDriverByName('GTiff')
ds = drv.Create('/Path/converted_numpy.tif', 5000, 6635,1,gdal.GDT_Float32)
outband=ds.GetRasterBand(1)

# Write array to raster
outband.WriteArray(array)
outband.SetNoDataValue(-999)
outband.FlushCache()

# Take care of spatial reference
ds.SetGeoTransform((414733,0.5,0,5953200,0,-0.5))
dsSRS=osr.SpatialReference()
dsSRS.ImportFromEPSG(25832)
ds.SetProjection(dsSRS.ExportToWkt())

Does anybody know if something's missing?
EDIT:
The Output of the query and reconverting is as follows:
cur.execute('''
  SET postgis.gdal_enabled_drivers TO 'GTiff';
  WITH temp as (
    SELECT ST_Union(rast) As file_rast
    FROM public.table) SELECT ST_AsGDALRaster(file_rast, 'GTiff') FROM temp;''')

for row in cur:
    np.rast = row[0].tobytes()
    with MemoryFile(np.rast).open() as dataset:
        data_array = dataset.read()
        print(data_array)
        max = np.max(data_array)
        print(max)
        print(type(max))
        print(type(max.item()))

Output: nan, <class 'numpy.float32'> and <class 'float'>


texts: Measure distance between vertices of polylines in QGIS for Desktop

texts:
I have this dataset of coordinates of electric lines in the form of vector polylines as shown in the image below, and I need to measure the distance between the first and last vertex of each line and all the other line vertices on the layer using PyQGIS. For example, take the first vertex of the first line and measure the distance between this one and all the others on the layer; then take the last vertex of the first line and measure all the distances, and repeat the same for the second line and so on.

I think that a good solution is to get the coordinates of the first and last vertices for each feature, save them as a point and then apply the method explained in this post: Measure distances between specific features in QGIS for Desktop to get the distances between points.


texts: Transform from WGS84 to EPSG 3785 (google maps)

texts:
I am a Java Developer and I was given this coordinates from a database

verbat_sistema_coordenadas  verbatim_SRS    datum   verbat_longitud verbat_latitud
UTM                         WGS84           WGS84   9947921          815887

And I tried the examples here https://stackoverflow.com/questions/7661/java-code-for-wgs84-to-google-map-position-and-back, but I get something like this when I try to use the MercatorTransform showed there

LatVerb           LonVerb               LatEPSG                        LonEPSG
815612.0         9948458.0           7.3267672551129115        66.40190209729147

I tried the formulas shown there and I get differnt result with each formula and none are close to what I need.  Then someone I am working with gave me an excel sheet with some formulas but the have a million variables and constants

Fi' Ni  a   A1  A2  J2  J4  J6  Alfa    Beta    Gamma   B(fi)   b   Zeta    Xi  Eta Sen h Xi    Delta Lambda    Tau

And they depend on the Timezone(meridian, maybe I am translating this wrong?), and the hemisphere, and they seem to work.  They gave me
-0.4544682953,-78.1736736245, which is the correct answer.

This is the excel sheet http://www.gabrielortiz.com/descargas/descarga.asp?Fichero=Conversion_UTM_Geograficas.xls

I tried using geotools, but it was way off too http://osgeo-org.1560.x6.nabble.com/convert-wgs84-coordinate-to-Google-Map-coordinate-system-EPSG-3785-with-geotools-and-java-td4983362.html


texts:
I have a graph with some data points:

From these data points I have created a table and calculated the Euclidean distance for each combination.

Now, I would like to calculate the semi variances for each combination. How do I proceed?

EDIT:

The expected results are shown in table below, however:


texts: Get Distance between two point on linestring using mysql Spatial

texts:
I have a linestring and two point on this linestrin according to my picture

i want distance between two point as base of linestring please hellp me


texts: How to perform enhancement operations on an satellite imagery in QGIS?

texts: I am a new user of QGIS please suggest how to perform image enhancement operations on an satellite imagery like ratioing, contrast stretching, edge engancement, filtering etc. in QGIS?


texts:
I would like for my points to cluster using L.markerClusterGroup.

I am receiving a 200 status; but nothing appears on my map.  I have commented out my markers which appear when uncommented.

http://jsfiddle.net/8mnyqx49/

/// var points = L.geoJson(geojsonfeatures, {
    ///filter: function(feature, layer) {
       /// return feature.properties.COUNTRY != "UNITED STATES";
  ///  },
   /// pointToLayer: function(feature, latLng) {
       ///  return L.marker(latLng)
         ///  .bindPopup(feature.properties.NAME);

var points = L.markerClusterGroup();

for (var i = 0; i < geojsonfeatures.length; i++) {
  var a = geojsonfeatures[i];
  var NAME = a[2];
  var marker = L.marker(new L.LatLng(a[0], a[1]), { NAME: NAME });
  marker.bindPopup(NAME);
  markers.addLayer(NAME);
}

map = L.map('myMap', {
  center: centerlatlng,
  zoom: 6,
  layers: [aLayerOne]
});


texts: ArcGIS - Only draw data in current layout view

texts:
I am using county data to create a map of a very localized area (about a city block). The data includes parcel, sewer, water, electric, stormwater, imagery, centerlines, etc.

Is there any way to force ArcMap to only draw the data that is within my current layout view? Each time I move the map an inch or add something it seems like I'm waiting a long time for it to draw everything out in the background.

I know I can select features manually and create a new layer or shapefile with them individually, however, there are quite a few dozen feature classes that need to be shown on this map.

This is something that would probably be useful in general and I'm sure it's a setting I have overlooked.


texts: How to interpolate integer rasters to floating point?

texts:
I have an reclassified raster, I need to convert it to floating point, the "float" tool does not suffice, because I need the output values to be spatially interpolated. For example, at the midpoint of zone 2 (halfway between borders of zone 1 and zone 3) I need the values to be 2.5, using the float tool these values would simply be 2.0

I am using ArcGIS 10.2


texts: Statistical tests for spatial line patterns?

texts:
There are plenty of tests for spatial point patterns that can be used to determine if points are distributed randomly or not, but are there any established tests for spatial line patterns? (I'm thinking straight lines, with just start and end point and no intermediate nodes.)

The data I want to analyze are OD (origin-destination) lines of human and animal movement. (Similar to the example in Clustering undirected lines.)

So far, one idea was to treat lines like 4D points and use point pattern tests but I'm uncertain if that is appropriate.

The ideal test would make it possible to determine whether there are clusters of lines or not.

Instinctively, I would say that many lines that start at the same origin but have all kinds of different destinations, should not be considered a cluster. On the other hand, many lines that run (close to) parallel for a longer time would be a cluster.


texts: OpenLayers.Format.GeoJSON throws an exception "SyntaxError: JSON.parse: unexpected character" for valid JSON

texts:
In my GeoExt client, I have the following line that retrieves JSON data from a web service into a FeatureStore:

var store = new GeoExt.data.FeatureStore({
        layer: route_layer,
        fields: [
            {name: "length"}
        ],
        proxy: new GeoExt.data.ProtocolProxy({
            protocol: new OpenLayers.Protocol.HTTP({
                url: "./php/pgrouting.php",
                format: new OpenLayers.Format.GeoJSON({
                    internalProjection: epsg_900913,
                    externalProjection: epsg_4326
                })
            })
        })
    });

When this code is executed, Firebug reports the error SyntaxError: JSON.parse: unexpected character, even though the JSON returned by the web service is valid (I've confirmed this against http://jsonlint.com). The response header has been configured accordingly with the line header('Content-type: application/json',true);.

Anyone who might have encountered this issue before, please assist. I've waded around, but no luck yet.

UPDATE: Sample JSON output is at http://pastebin.com/LPLjBDgp. It's valid JSON, but the error SyntaxError: JSON.parse: unexpected character persists. I've even updated my html page by adding the following line <meta http-equiv="Content-Type" content="text/html; charset=UTF-8"/> in an attempt to restrict character encoding.


texts: Rotating image by specific angle in QGIS plugin

texts:
I am using the Freehand Raster Georeferencing code in a widget to help size and locate site plans on aerial GE backdrops. I need to find a way for a user to input a specific angle of rotation and the site plan snaps to that angle. Below is a screenshot of my current set-up with comments.

In addition, is there a way to have a Qt-widget that dynamically changes the rotation angle when the user manually rotates the image. (Note: this is not my main question)

Below is the georeferencing rotation class:

class RotateRasterMapTool(QgsMapToolEmitPoint):
def __init__(self, iface):
    self.iface = iface
    self.canvas = iface.mapCanvas()
    super(RotateRasterMapTool, self).__init__(self.canvas)

    self.rasterShadow = RasterShadowMapCanvasItem(self.canvas)

    self.rubberBandExtent = QgsRubberBand(self.canvas, QgsWkbTypes.LineGeometry)
    self.rubberBandExtent.setColor(Qt.red)
    self.rubberBandExtent.setWidth(1)

    self.reset()

def setLayer(self, layer):
    self.layer = layer

def reset(self):
    self.startPoint = self.endPoint = None
    self.isEmittingPoint = False
    self.rubberBandExtent.reset(QgsWkbTypes.LineGeometry)
    self.rasterShadow.reset()
    self.layer = None

def canvasPressEvent(self, e):
    self.startY = e.pos().y()
    self.endY = self.startY
    self.isEmittingPoint = True
    self.height = self.canvas.height()

    self.isLayerVisible = isLayerVisible(self.iface, self.layer)
    setLayerVisible(self.iface, self.layer, False)

    rotation = self.computeRotation()
    self.showRotation(rotation)

def canvasReleaseEvent(self, e):
    self.isEmittingPoint = False

    self.rubberBandExtent.reset(QgsWkbTypes.LineGeometry)
    self.rasterShadow.reset()

    rotation = self.computeRotation()
    self.layer.setRotation(self.layer.rotation + rotation)

    setLayerVisible(self.iface, self.layer, self.isLayerVisible)
    self.layer.repaint()

    self.layer.commitTransformParameters()

def canvasMoveEvent(self, e):
    if not self.isEmittingPoint:
        return

    self.endY = e.pos().y()
    rotation = self.computeRotation()
    self.showRotation(rotation)

def computeRotation(self):
    dY = self.endY - self.startY
    return 90.0 * dY / self.height

def showRotation(self, rotation):
    center, originalRotation, xScale, yScale = self.layer.transformParameters()
    newRotation = rotation + originalRotation
    cornerPoints = self.layer.transformedCornerCoordinates(center, newRotation, xScale, yScale)

    self.rubberBandExtent.reset(QgsWkbTypes.LineGeometry)
    for point in cornerPoints:
        self.rubberBandExtent.addPoint(point, False)
    # for closing
    self.rubberBandExtent.addPoint(cornerPoints[0], True)
    self.rubberBandExtent.show()

    self.rasterShadow.reset(self.layer)
    self.rasterShadow.setDeltaRotation(rotation, True)
    self.rasterShadow.show()

Any ideas as to how to integrate a user-specific rotation angle input with this dock.

Any ideas as to how I could take the rotation angle and show it on a Qt-widget?


texts: QGIS3 shows coordinate in geographic coordinates though QGIS project and layer are both in projected CRS

texts:
The CRS of the data I am using is in geographic CRS (WGS 84 EPSG:4326), but I want to work in projected CRS which for my case is in UTM 32N, EPSG 25832. I am working on modifying polygons with the advanced digitizing panel which only works on projected CRS.

I am able to project the layer and also set the project CRS to EPSG 25832, but the coordinates (marked in screen-shot) is always shown in lat/lon. Also when I use the digitizing toolbox to modify the position of the feature, the unit shown in the panel are also lat/lon, which is strange because the CAD tool does not support geographic coordinates (I cant use the tool when the project CRS is set to geographic CRS).

Any idea what the issue is? I am working with QGIS 3.4.1, and in QGIS 3.2 I have the same problem. I tried QGIS 2.18 where I don't have this issue.


texts: Why can't I multiply values from column named 'length' using field calculator via PyQGIS?

texts:
I want part of my plugin to multiple value of length and add it to table.Basing on this question: How to use QGIS field calculator in python? , I made script that works fine for some filed: I can multiple values of another field (integer) but I can't multiple values ot the 'length' field (double type). I don't know why, but when I'm trying to do : 'length *2' the result is NULL. The 'length' is calculated also inside the plugin.

I'm not sure if the question is understandable, hope code snippets will help:

This is how I calculate length to the 'length' field:

layer = QgsMapLayerRegistry.instance().mapLayersByName("memory:temp_layer") [0]

layer.dataProvider().addAttributes([QgsField("length", QVariant.Double,len = 6, prec = 2)])
layer.updateFields()

idx = layer.fieldNameIndex('length')

layer.startEditing()

for feat in layer.getFeatures():
    lengths = feat.geometry().length()
    layer.changeAttributeValue(feat.id(), idx, lengths)

layer.commitChanges()

This is how I want to multiple:

layer = QgsMapLayerRegistry.instance().mapLayersByName("memory:temp_layer")[0]

layer.startEditing()

layer.dataProvider().addAttributes([QgsField('n', QVariant.Double,)])
layer.updateFields()

idx = layer.fieldNameIndex('n')

e = QgsExpression ("length * 2")
e.prepare(layer.pendingFields())

for f in layer.getFeatures():
    f[idx] = e.evaluate(f)
    layer.changeAttributeValue(f.id(),idx, e.evaluate(f))

layer.commitChanges()

I came up with idea that the problem may be that lengths type is float while Qvariant type of field is double. Even if, I don't know how to handle this.

I can multiple 'length' field using standard field calculator, not via Python script and the result is fine.

Any ideas?


texts: Displaying attribute tables with merged cells in ArcGIS Desktop?

texts:
I am new to ArcGIS Desktop 10.1. I would like to display a table providing information about the Region, Province and district in the layout view of ArcMap.

For example, this is my attribute table:

I want to display the following table in the layout view of ArcMap:


texts:
I am currently looking for a DEM containing some drumlins to perform some morphometry on it. It should be a free one, or at least to be downloaded free of charge.

I have already tried to find a free DEM for Drumlin spots in the Alps foreland in southern Germany, which failed. The point is that I need a DEM with a higher resolution than 30 meters/pixel (which would equal the 1 arcsec of SRTM). Minimum would be 10 meters, better 5 or even 1 meter per pixel. Do you know any options, at least containing drumlins?


texts: Getting MemoryError making hexagonal lattice over polygon in ArcPy?

texts:
I am creating a hexagonal lattice that I will then use to create Thiessen Polygons. We are required to make a lattice first. I can't figure out how to get rid of the MemoryError that looks like this:

Traceback (most recent call last):
  File "C:\Users\shehn\Documents\ArcGIS\PyScripts\lab6ed.py", line 31, in
    sh.append(arcpy.PointGeometry(p))
MemoryError

Failed to execute (Hex1).

Here is the code I'm using

    # Import arcpy module
import arcpy

# Local variables:
lyr = arcpy.GetParameterAsText(0)
b = arcpy.mapping.Layer(lyr)
file = arcpy.GetParameterAsText(1)
spacing = int(arcpy.GetParameterAsText(2))
if(not file.endswith(".shp")) : file += ".shp"
e = b.getExtent()

xamin = e.XMin
yamin = e.YMin
xamax = e.XMax
yamax = e.YMax
deltaX = xamax - xamin
deltaY = yamax - yamin
half = spacing / 2

numpx = int(deltaX / spacing) + 2
numpy = int(deltaY / spacing) + 2

p = arcpy.Point() # Create a "Point" object
sh = []

for j in range(numpy):
    p.Y = yamin + spacing * j
    for i in range(numpx):
        p.X = xamin + spacing * i
        if (j % 4 == 1 or j % 4 == 2) : p.X += half
        sh.append(arcpy.PointGeometry(p))

temp = "C:\\Temp\\t.shp"
arcpy.CopyFeatures_management(sh, temp) # Create a temp shape file

arcpy.Clip_analysis(temp, lyr, file, "")
arcpy.Delete_management(temp)
# Create a layer from the shape file and add it to the map
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layer = arcpy.mapping.Layer(file)
arcpy.mapping.AddLayer(df, layer)

# Delete unused objects
del layer
del df
del mxd


texts:
I am new to ArcGIS and stuck with "Error HRESULT E_FAIL has been returned from a call to a COM component." error while converting gdb to geodatabase through Arcgis .net SDK verstion (Arcgis Desktop version 10.3.1). Following are my complete code

using System;
using System.Collections.Generic;
using System.Text;
using ESRI.ArcGIS.esriSystem;
using System.Configuration;
using ESRI.ArcGIS.DataManagementTools;
using ESRI.ArcGIS.Geoprocessor;
using System.Runtime.InteropServices;

namespace Drainage_geodbCreator
{
    class Program
    {

        private static LicenseInitializer m_AOLicenseInitializer = new Drainage_geodbCreator.LicenseInitializer();

        [STAThread()]
        static void Main(string[] args)
        {
            String MXDpath = ConfigurationManager.AppSettings["MXDpath"];
            String geoDBPath = ConfigurationManager.AppSettings["geodatabsePath"];
            int mycheck = Convert.ToInt32(ConfigurationManager.AppSettings["geodatabaseversion"]);
            ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
            //ESRI License Initializer generated code.
            m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeStandard },
            new esriLicenseExtensionCode[] { });

            CreateRuntimeContent _CreateRuntimeContent = new CreateRuntimeContent();
            _CreateRuntimeContent.in_map = MXDpath;
            _CreateRuntimeContent.output_folder = geoDBPath;
            _CreateRuntimeContent.extent = "DEFAULT";
            _CreateRuntimeContent.options = "FEATURE_AND_TABULAR_DATA";
            _CreateRuntimeContent.optimize = "NON_OPTIMIZE_SIZE";

            bool rst = false;
            try
            {
                ExceuteGeoProcesser(_CreateRuntimeContent, out rst);
            }
            catch (COMException ce)
            {

            }

            //ESRI License Initializer generated code.
            //Do not make any call to ArcObjects after ShutDownApplication()
            m_AOLicenseInitializer.ShutdownApplication();
        }

        static Boolean ExceuteGeoProcesser(IGPProcess _GPProcess, out Boolean ProcessResult)
        {

            Boolean Rslt = false;

            try
            {
                Geoprocessor _Geoprocessor = new Geoprocessor();
                Guid g = new Guid("5374EC4C-1AA2-4829-A811-DE624ECEC23F");
                Type t = Type.GetTypeFromCLSID(g);
                _Geoprocessor = Activator.CreateInstance(t) as Geoprocessor;
                _Geoprocessor = new Geoprocessor();

                _Geoprocessor.Execute(_GPProcess, null);

                ProcessResult = Rslt;

                return true;

            }
            catch (Exception ex)
            {
                ProcessResult = Rslt;
                return false;
            }

        }
    }
}


texts: LASToMultiPoint_3d creating empty files?

texts:
Trying to use the arcgis python command
LAStoMultipoint_3d as seen here.

I'm running the code snippet against one of my LAS files.  It creates a .shp, sbx, shp.xml, shx and dbf file.  However, the files are empty.

The code is as follows:

#import script modules
import arcpy
from arcpy import env

arcpy.CheckOutExtension("3D")
env.workspace = "D:\myFolder\"

arcpy.LASToMultipoint_3d("D:/myFolder/myFile.las",
                         "D:/myFolder/Output" , 1.5)

print 'done'

The code executes and creates empty files.  I have tried this with a half-dozen different files varying in size and complexity.


texts:
I have a custom made plugin for QGIS, which uses Java JDK (running on Win7x64). If the JAR file is run through command prompt, everything works great, but when the same command is used in plugin, it shows this error:

Exception in thread "main"
java.lang.UnsatisfiedLinkError:
com.sun.imageio.plugins.jpeg.JPEGImageReader.initReaderIDs(Ljava/lang/Class;Ljava/lang/Class;Ljava/lang/Class;)V
at com.sun.imageio.plugins.jpeg.JPEGImageReader.initReaderIDs(Native Method)
at com.sun.imageio.plugins.jpeg.JPEGImageReader.<clinit>(JPEGImageReader.java:91)
at
com.sun.imageio.plugins.jpeg.JPEGImageReaderSpi.createReaderInstance(JPEGImageReaderSpi.java:85)
at javax.imageio.spi.ImageReaderSpi.createReaderInstance(ImageReaderSpi.java:320)
at javax.imageio.ImageIO$ImageReaderIterator.next(ImageIO.java:529)
at javax.imageio.ImageIO$ImageReaderIterator.next(ImageIO.java:513)
at com.geokods.App.main(App.java:262)

I fount that QGIS uses different Java version:

java version "1.7.0_67"
Java(TM) SE Runtime Environment (build 1.7.0_67-b01)
Java HotSpot(TM) 64-Bit Server VM (build 24.65-b04, mixed mode)

vs. CMD:

java version "1.8.0_45"
Java(TM) SE Runtime Environment (build 1.8.0_45-b14)
Java HotSpot(TM) 64-Bit Server VM (build 25.45-b02, mixed mode)

How can I fix this?

UPDATE:

Problem is partially solved by appending command with:

-version:1.8.0_45

Now plugin works as it should. But still, there should be a way to overwrite the old Java version.

UPDATE 2:

Tried to set JAVA_HOME by adding set JAVA_HOME=C:\Program Files\Java\jdk1.8.0_45 to bin\qgis.bat. Did not help.


texts: Copy/Paste a row (duplicate an object) in the same attribute table

texts: When editing in QGIS in some shape attribute tables I can copy/paste a row with data, but in some shapes when I do this it adds NULL data rows, I know there is a solution but I failed to find out or to find here. Please help...


texts: How to combine date and time in field calculator in arcgis

texts:
I want to simulate the travel path based on data points generated by my farm tractor, which is equipped a GPS receiver (I want to know where it starts and where it ends). Unfortunately the generated data points are not labelled sequential IDs.

In ArcGIS, I have two columns Date and Time on the attribute table that I can refer to as sequence and then use points to lines to generate lines.

Both Date and Time are string field:

Maybe I should combine two columns into one? I've tried few times inputting Time as Line Field and Date as Sort Field but turns out warning001059.


texts: Turning a JSON data to a GeoJSON data. For Leaflet


texts: How to debug PostGIS connection issue from MapServer?

texts:
I have a PostGIS connection string that works fine in GDAL, as I get a raster response when running

gdalinfo "PG:host=... port=5432 user='...' password='...' dbname='weather_models' schema='rasters' table='wpc_snow_50pct_1551268800' column='rast' mode=1 where='timestamp=\'2019-03-02 00:00:00\''"

On this very same server, I have cgi-mapserver running, and the output from /usr/lib/cgi-bin/mapserv/ -v is

MapServer version 7.2.1 OUTPUT=PNG OUTPUT=JPEG OUTPUT=KML SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=SVG_SYMBOLS SUPPORTS=RSVG SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER SUPPORTS=SOS_SERVER SUPPORTS=FASTCGI SUPPORTS=THREADS SUPPORTS=GEOS SUPPORTS=PBF INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

The problem is, when I use that same PG:... connection string as the DATA for a LAYER in MapServer, I get an unhelpful error:
Unable to access file. Corrupt, empty or missing file '<snip of PG connection string>' Couldn't establish a database connection.

Despite turning DEBUG level to 5, and CPL_DEBUG to ON in the MAP block of my mapfile, I don't get any helpful debug information whatsoever in my error log -- just the connection string and "Couldn't establish a database connection."  I have tried a variety of where clauses and other tables and get the same result.  I have no idea if the actual error is coming from GDAL or MapServer, despite all those strings working fine when running them directly via GDAL.

Here is the mapfile:

MAP
 CONFIG "MS_ERRORFILE" "/map/error.txt"
 CONFIG "CPL_DEBUG" "ON"
 DEBUG 5
 PROJECTION
  "init=epsg:4326"
 END
 LAYER
  NAME wpc_snow_50pct
  TYPE RASTER
  STATUS ON
  VALIDATION
   "timestamp" ".*"
  END
  PROJECTION
   "init=epsg:4326"
  END
  METADATA
   "wms_title" "wpc_snow_50pct"
   "wms_srs"   "EPSG:4326"
   "wms_extent" "-108.5 41 -104.5 36"
   "wms_enable_request" "*"
  END
  DATA "PG:host=... port=5432 dbname='weather_models' user='...' password='...' schema='rasters' table='wpc_snow_50pct_1551268800' column='rast' mode=1 where='timestamp=\'2019-03-02 00:00:00\''"
  PROCESSING "BANDS=1"
 END
END


texts: What is Stopping Geometric Network Trace While No Barrier Specified

texts:
I have an Electric geometric network which has a Generator as a Source and Flow has been set up like below image:

What is stopping the tracing when no Barrier been specified? For example:

As you can see there are no Barriers specified here, so why is the Trace not continuing until the Source or last sink at Upstream or Downstream Trace?


texts: Configuring Geoserver to use a specific postgis geometry column as the default

texts:
I have a PostGIS table with two geometry columns, both multipolygon but one with SRID 4326, the other 3857. When publishing the table in Geoserver, I realize that Geoserver defaults to using the first geometry column in the list of columns, but I was wondering if there is a way to configure Geoserver to use a specific geometry column?

I’m not talking about configuring the SLD to render using a specific geometry column, ie:

<se:PolygonSymbolizer>
<se:Geometry>
<ogc:PropertyName>geom3857</ogc:PropertyName>
</se:Geometry>
…
</se:PolygonSymbolizer>

I’m talking about having Geoserver use a specific geometry column as the native SRID and default projection.


texts: Spatially joining polygon fields onto point layer in QGIS?

texts:
I have gridded depth points in CSV format that I'm trying to add other map data to, in additional columns. I've tried joining fields, hoping that I could link the Shapefile data to the csv, but when I try to display the csv by graduated symbols based on the new fields, it goes blank. I used Nathan's tips from here, but imported the csv as csv since I want to add map info to that, rather than adding csv info to a map.

See this picture below for the sea bottom temperature & salinity voronois by month, underneath the csv points (a sample subset), with the joins already made. I added columns on the csv for each of the fields i'm intending to append.

I also tried spatial joins as per this however it says the layers aren't the same CRS (I checked, they are), then gets to 15% and hangs.

I also also tried spatial join through mmqgis but it hit a python error. I can try this again and look to debug the problem if this is the right way to go.

In short: appending polygon data as fields on a points-grid csv: what's the best way in QGIS 2.0?


texts: Determining coordinate system of image to be georeferenced?

texts:
Possible Duplicate:
  .shp with unknown coordinate system

I need to capture spatial information off of an old map and I can't determine the coordinate system of the image.  I also need to examine/cross-reference some information via overlay.  The subject area is Kirkland Lake Ontario, the scale of the map is 1 inch to 400 feet, and an example of the coordinates shown is 86,000E - 90,000N (on the southwest corner of the map).

I took a few guesses at the coordinate system (as did my co-worker) and we haven't had any success.  We figured it might be something of the old NAD27 variety.  I also tried using the road, rail, utility line and river layers from the province to georeference the old map, but I'm not having good results that way either.  Accuracy is preferred too, so using the grid to georeference this map would be the preferred method.  I fear this map may have been generated with some sort of custom coordinate sytem that the local mining companies were using many decades ago.

Any ideas?


texts: Using CJK character orientation in label expression of ArcMap?

texts:
I'd like to set label expressions by VBscript in ArcGIS, here comes the problem.
Since labels are constructed by Chinese and English, so I set the label expression like:

If CJK character is checked, Chinese labels would ideally arranged whereas English labels would not.

(1) CJK character is not checked:

(2) CJK character is checked:

I want the labels to be arranged like Chinese in (2) and English in (1).


texts: Is it valid to adjust quantile buckets to "nice" numbers for better readability?

texts:
I'm making a map about terrorist attacks and fatalities in CartoDB. For data clustering I'm using quantiles. The calculated buckets are 16, 32, 50, 97, 700. Is it okay to adjust the quantiles to make them more readable, like 15, 30, 60, 90, 700?

Personally, I think it doesn't matter, since the data clustering is arbitrary and the map would look totally different, if I had used a different clustering method.


texts: Calculating shortest distance between polygons?

texts:
I have many polygons stored in a table. I want to filter out polygons where they separate too far away (say 1km).

My approach is:
calculate the minimum distance between a particular polygon to other polygons, and store the value in a new column in the table.

Then filter out polygon with minimum distance larger than 1 km.

I have tried to use:

ALTER TABLE "int_20190124" ADD "nearestDistance" float8;
INSERT INTO "int_20190124"("nearestDistance")
VALUES
(
SELECT ST_Distance("int_20190124".geom, "int_20190124".geom)
FROM "int_20190124"
ORDER BY ST_Distance("int_20190124".geom, "int_20190124".geom)
LIMIT 1
)

I am very new to postgis.


texts: Weird coordinates after converting numpy array to raster

texts:
I have the NumPy array and its bounding box coordinates. I have tried to convert it into raster using rasterio, based on this answer, and it did save it as raster, but when I use rasterio.show the coordinates are very wrong.
This is the script I have used:
bbox_coords_wgs84=[-101.7359960059834, 20.21904081937658, -100.5717967351885, 20.8312118894487]

#variables for the projection:
minx=bbox_coords_wgs84[0]
maxy=bbox_coords_wgs84[3]
pixel_size= 10

#according to the post on GIS SO:

import rasterio
from rasterio.transform import from_origin

transform=from_origin(minx,maxy,pixel_size,pixel_size)
crs_img='EPSG:4326'

with rasterio.open('test1.tif',
                    'w',
                    driver='GTiff',
                    height=ndvi.shape[0],
                    width=ndvi.shape[1],
                    count=1,
                    dtype=ndvi.dtype,
                    crs=crs_img,
                    nodata=None, # change if data has nodata value
                    transform=transform) as dst:
        dst.write(ndvi, 1)

#display the results:

from matplotlib import pyplot
from rasterio.plot import show

src = rasterio.open('test1.tif')
show(src)

As you can see, the numbers are absolutely not the correct coordinates.
My end goal: to be able to reproject the NumPy array into WGS84 correctly.
*This post relates also to this post


texts: Any ZXY satellite basemaps covering New Zealand at zoom 19?

texts:
I need a good satellite base map to use with CartoDB.  It has to cover all of New Zealand at zoom 19, [edit:] and be as recent as possible so we can show recently-built properties.  It doesn't have to be free. (In fact, we're rather expecting to have to pay for it.)

I'm using ArcGIS's http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x} at the moment, and while it's good at zoom 19 in some places...

... in other places it only goes to 18, and it can get pretty grainy...


texts: How to merge multiple file geodatabases into one file geodatabase?

texts:
I have created a file geodatabase that has several feature data sets, feature classes, and domains.  I have created a map with the file geodatabase added to it, converted the map to ARCPAD AXF format and copied the map to several Junos changing the name for each Juno unit.  Data is being successfully collected.

My question is: How do I merge all these seperate geodatabases from the Junos into the central file geodatabase that I created originally?


texts: ntv2 datum transformations-problem on file output

texts:
I am working with qgis 2.18.27
I have installed ntv2 datum transformations and it worked fine, but yesterday it didn't work and i am not finding the way to resolve the issue (explained below)

I am using the portugal ntv2 grids (direct and inverse vector transformation)

The plugin after processing ouputs a blank shapefile.

Can anyone help me?


texts: simplify tiled geometries while preserving anything on the edge

texts:
I have a dataset of polygons split into evenly sized tiles. How can I go about simplifying them, but not along the tile edges, so that there aren't any gaps between continuous features when the tiles are aligned with one another? For instance, here are two unsimplified tiles:

After simplification, a "seam" forms between the two halves of Antarctica:

Here are Gists of the non-simplified and simplified GeoJSON files. I'm currently splitting the original geometries into tiles and then simplifying, and one solution appears to be reversing that order (simplifying and tiling afterwards), but that's unfortunately not feasible for architectural reasons. A PostGIS solution is preferred, but alternatives are welcome!


texts: Convert folder of XYZ png tiles to TMS Tiles

texts:
I have elevation/terrain tiles downloaded from https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png  I need to convert them to TMS Tiles.

Is there a python script that can do the re-numbering?


texts: How to use Python implemented GDAL utilities via OSGeo4W Shell?

texts:
Following installation of QGIS v 3.10.1 on a W10 machine, I am automatically hooked up with OSGeo4W Shell. From there, I can directly invoke gdal utilities such as gdalinfo and gdal_translate. However, other utilities fail me, such as gdal_edit and pct2rgb, even though they are listed alongside each other at the gdal raster programs website.
It appears that gdal_edit and the likes are implemented in Python, quoting from OSGeo4W site:

To get the commandline utilities implemented in Python (like gdal_merge, and rgb2pct), install gdal-python.

Reading up on gdal-python says:

The standard scripts, like rgb2pct.py, are placed in osgeo4w\bin and include .bat wrapper scripts

This is as far as I get. I can't find the mentioned location, nor any .bat wrapper scripts. The nearest I can get is: c:\Program Files\QGIS 3.10\apps\Python37\Scripts\gdal_edit.py but I am lost how to employ this script to my needs.
How can I use gdal_edit?


texts: Connecting from ArcGIS Desktop to virtual machine's ArcGIS Server?

texts:
I have made an OpenVPN connection to an "inner network", my school's private network. I can use the remote desktop application after the connection, so it's surely work. The ArcGIS Server is in a virtual machine, in a Windows Server 2008 R2.

The Server's version is 10.1 (I did not have refreshed that yet), and the ArcGIS Desktop version is 10.2. The Desktop application is installed on my computer. I have added myself as administrator on the ArcGIS Server Manager.

How should I connect to the server (i.e. add that GIS server to the ArcCatalog)?

I have written in the URL:

http://nameofthecomputerwithArcGISServer:6080/arcgis

But it says it cannot connect to the server.


texts: How to pan-sharpen Landsat images for image-classification in GRASS?

texts:
I'd like to learn how to do an unsupervised classification of a Landsat scene using i.cluster > i.maxlik in GRASS using pan-sharpened 15m-resolution images (the example given in its Wiki uses the ordinary 30m-resolutions).

I tried using i.pansharpen to create the pan-sharpened images first, but i.pansharpen can only output 3 channels that can be combined using d.rgb or r.composite. As far as I know, the image-classification process requires complete multispectral bands 1 to 7. How can I produce separate bands (band 1 to band 7), pan-sharpened using its 15m-resolution band 8 image, before inputting them to the image-classification process?

I've found a paper that demonstrates how they did this; basically they used some kind of Principal Component Analysis to somehow merge the 30m-multispectral bands with the 15m-panchromatic band. The exact quote would be:

  "The method first resamples the 30-meter multispectral image to 15 meters.  It then computes all six principal components of the multispectral image  (We deleted the thermal band due to is courser resolution).  Then the histogram of the panchromatic band (15-meter resolution) is rescaled to match the histogram
  of the first principal component of the 30-meter image and the first component is replaced with the rescaled panchromatic band.  This is justified because the first principal component represents overall brightness in a way that is similar to the broad spectral band of the panchromatic image.  After replacement, then the six components are transformed back into the original data space, rendering an improvement in spatial resolution."

First off, the paper showed no algorithm/formulas whatsoever. I have no idea how to turn the above quote into a corresponding mathematical formula. I realized that I can use i.pansharpen with the PCA algorithm instead of the usual Brovey or IHS - but still - the output will only be 3 channels of red, blue and green - which I sadly have no clue on how to use them for image-classification..

So, before I even attempt to crack my head trying to write a new PCA algorithm manually, can anybody help me to point to an easier and better way to perform image-classification on pan-sharpened Landsat images? I mean - there should be an easier way, right? I feel like I'm missing something simple.

If the only route left is to write my own script, can you guys please point me to something out there that resembles an example of what I'm trying to do?

Any help is very much appreciated!


texts: QGIS server + spawn-fcgi : where can I get informations when an error occurs?

texts:
I am new to GIS, QGIS and spawn-fcgi.
Thus, the answer to my question may be obvious.... but I couldn't figure it out (or express it in a way an accurate way such as I find accurate answers!)

Here it is :
I retrieved a docker Image of QGIS server (3.0) that is served through a (also docker-hosted) nginx server and spawn-fcgi. The configuration looks very standard, and exactly follows QGIS Getting Started tuto.

The .qgs project provided as as example runs well : I perform some WMS requests on it, and retrieve the correct data.

But when I switch to one of my projects (whose WMS capabilities have also been set), I get an error.
Whatever the request is, it keeps failing:

http://172.....:8080/qgis?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&Layers=mylayer&BBOX=1,2,1.1,2.1&...
http://172.....:8080/qgis?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities
http://172.....:8080/qgis?SERVICE=WMS
http://172.....:8080/qgis?SERVICE=toto

=> 500

<ServerException>Project file error</ServerException>

I followed a few tips provided on this post , but though I fixed some volume mounting problem, I still have my pb.

The project is present under the expected directory (specified by QGIS_PROJECT_FILE):

$ qgis exec -it mycontainer bash
qgis@mycontainer $ ls ${QGIS_PROJECT_FILE}
/data/myproject.qgs
qgis@mycontainer $

Also, when I consult logs (with debug activated and log level set to all :'0'), I see the WMS requests :

[21][09:12:17] ******************** New request ***************
[21][09:12:17] REMOTE_ADDR: 172.18.0.1
[21][09:12:17] CONTENT_TYPE:
[21][09:12:17] HTTP_USER_AGENT: Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:54.0) Gecko/20100101 Firefox/54.0
[21][09:12:17] REQUEST:getCapabilities
[21][09:12:17] SERVICE:1.3.0
[21][09:12:17] Request finished in 8 ms

The ngninx logs show the 500 error... and that's all.

And here's my point : how do you get informations about qgis server errors?

Looks like I should have a look to spawn-fcgi logs... but I am quite lost here : I eventually have no clue on where I should look for them, or whatever informations about how to understand (and fix!) my problem...
So any tips are welcome!


texts: Exclude geographic extent of an OSM file

texts:
Osmium and Osmosis command tools have the capability of extracting a geographic extent from an OSM file.

osmosis
--read-xml test.osm --tee 1
--bounding-box left=15 top=46 --write-xml extract_test.osm

osmium extract -p city.geojson test.pbf -o extract_test.pbf

What I need is to extract region from an OSM file like a city and after some edits, merge it with the OSM file with the hole (excluded the city region). How can I exclude or clip a geographic extent from original file?


texts:
How do I use zonal statistics in the python console to get CSV/dbf output?

I have used the code like this:

# Python code for Zonal statistics
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
from qgis.analysis import QgsZonalStatistics
#specify polygon shapefile vector
polygonLayer = QgsVectorLayer(' E:/Trash/BhagaP.shp', 'zonepolygons', "ogr")
# specify raster filename
rasterFilePath = 'E:/Trash/Test48P.img'
# usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1)
zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, 'pre-', 1, QgsZonalStatistics.Mean)
zoneStat.calculateStatistics(None)


texts: ArcObjects error storing long integer value to Oracle

texts:
I am trying to save a long integer value to an Oracle field. I have created a field of long integer type using ArcCatalog in a standalone table. The below code breaks at row.set_Value(idx, val) . Error message is:

  "The value type is incompatible with the field type.[EMP_ID]"

It looked such a trivial problem but now I am stuck. Any ideas how to resolve this?

IMxDocument doc = ArcMap.Application.Document as IMxDocument;
IMap map = doc.FocusMap;
IFeatureLayer layer = map.get_Layer(0) as IFeatureLayer;
IFeatureClass fc = layer.FeatureClass;
IWorkspace ws = (fc as IDataset).Workspace;

IFeatureWorkspace fws = ws as IFeatureWorkspace;
ITable tbl = fws.OpenTable("TestTable");
int idx = tbl.FindField("EMP_ID");
IRow row = tbl.CreateRow();
long val = 123412341234;
row.set_Value(idx, val);
row.Store();


texts:
I have a geo-referenced pdf that I'm transforming into .mbtiles or .tiff using

gdal_translate 00989I11.PDF 00989I11.tiff

The problem is that the result does not have a white background as I would expect. I'm a bit new to all the gdal manipulation but my question is how could I change all the transparent pixels into white pixels.

PDF looks like this. source here

The resulting .tiff looks like this source here:

The gdalinfo of the .pdf is:

gdalinfo 00989I11.PDF
Driver: PDF/Geospatial PDF
Files: 00989I11.PDF
Size is 806, 1237
Coordinate System is:
PROJCRS["",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.2572221,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6269]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",35.2372777777778,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-120.642611111111,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",45,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",33,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["Inch",0.0254000508001],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["Inch",0.0254000508001],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["Inch",0.0254000508001]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["Inch",0.0254000508001]]]
Data axis to CRS axis mapping: 1,2
Origin = (-1807172.000000000000000,2042985.000000000000000)
Pixel Size = (3334.000000000000000,-3335.000000000000000)
Metadata:
  AUTHOR=
  CREATION_DATE=D:20190718071421-04'00'
  CREATOR=MicroStation 8.11.9.867 by Bentley Systems, Incorporated
  NEATLINE=POLYGON ((-1487388.02470535 -1654179.3669621,563048.661741154 -1654193.12139478,563073.777508997 1615063.22966106,-1487370.34966788 1615075.17034221,-1487388.02470535 -1654179.3669621))
  PRODUCER=Adobe PDF Library 9.0
  TITLE=00989I11
Corner Coordinates:
Upper Left  (-1807172.000, 2042985.000) (121d 9' 6.29"W, 35d42'19.70"N)
Lower Left  (-1807172.000,-2082410.000) (121d 8'43.15"W, 34d45'28.56"N)
Upper Right (  880032.000, 2042985.000) (120d23'40.84"W, 35d42'22.90"N)
Lower Right (  880032.000,-2082410.000) (120d23'52.11"W, 34d45'31.72"N)
Center      ( -463570.000,  -19712.500) (120d46'20.58"W, 35d13'57.63"N)
Band 1 Block=806x1 Type=Byte, ColorInterp=Red
  Overviews: 403x619, 202x310, 101x155
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 403x619, 202x310, 101x155
Band 2 Block=806x1 Type=Byte, ColorInterp=Green
  Overviews: 403x619, 202x310, 101x155
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 403x619, 202x310, 101x155
Band 3 Block=806x1 Type=Byte, ColorInterp=Blue
  Overviews: 403x619, 202x310, 101x155
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 403x619, 202x310, 101x155
Band 4 Block=806x1 Type=Byte, ColorInterp=Alpha
  Overviews: 403x619, 202x310, 101x155


texts: Reading data from nearest point and joining it to another point layer?

texts:
I have a set of points (villages) and another set of grid points that have a distance value. I need to identify the nearest grid point to the village points and then read the distance value and join it to the village.

Can you please tell me how to do this with basic license, and if not, with Standard license?


texts: Uncheck a layer (visibility) on the legend panel using PyQGIS?

texts: I'm trying to uncheck a layer on the layer panel using PyQGIS such that it is not visible on the project without removing it (like with removeMapLayer).


texts: Simplifying a long if condition statement and adding status report?

texts:
def Assign(csdname, cdname, ccsname, zone):
if csdname == "Oshawa" or csdname == "Whitby":
    return "70"
elif csdname == "Toronto":
    return "71a"
elif csdname == "Mississauga" or csdname == "Brampton":
    return "71b"
elif csdname == "Hamilton" and cdname == "Hamilton":
    return "72"
elif csdname == "Ottawa" or csdname == "Vanier" or csdname == "Rockcliffe Park":
    return "73a"
elif csdname == "Kanata" or csdname == "Nepean" or csdname == "Cumberland" or csdname == "Gloucester":
    return "73b"
elif csdname == "Windsor" or csdname == "Tecumseh" or csdname == "St. Clair Beach" or csdname == "Sandwhich West":
    return "74a"
elif csdname == "London":
    return "74b"
elif csdname == "Niagara Falls":
    return "75a"
elif csdname == "St. Catharines":
    return "75b"
elif csdname == "Burlington" or csdname == "Oakville":
    return "76a"
elif csdname == "Kitchener" or csdname == "Cambridge" or csdname == "Waterloo":
    return "76b"
elif csdname == "Sudbury":
    return "77a"
elif csdname == "Sault Ste. Marie":
    return "77b"
elif zone == "78a":
    #78 Northern Ontario A
    return "78a"
elif csdname == "Sudbury" or csdname == "Timiskaming":
    #78 Northern Ontario B
    return "78b"
elif cdname == "Algoma" or cdname == "Cochrane" or cdname == "Kenora":
    #78 Northern Ontario C
    return "78c"
elif csdname == "Thunder Bay":
    return "79b"
elif zone == "West 85th Meridian":
    return "79a"
else:
    return "70"

Assign(!CSDNAME!, !CDNAME!, !CCSNAME!, !ZONE!)

Trying to simplify the Python script above but not sure how to best do it. Ideally, I'd also like to figure out how to print which of the if conditions are not met for manual review later.


texts: Reading and plotting shapefiles with a custom function in R?

texts:
I have this piece of code that want to turn into a function

r <- readShapeSpatial("roads")
b <- readShapeSpatial("buildings")
lc <- readShapeSpatial("landcover.shp")

here is the preparation

laer <- function(x,y,z) {
  r <- readShapeSpatial(x)
  b <- readShapeSpatial(y)
  lc <- readShapeSpatial(z)
  plot(lc)
  plot(b,add=T)
  plot(r,add=T)
}

Then i use the function that i just created

laer("railways.shp","buildings.shp","natural.shp ")

and i get this

Error in getinfo.shape(filen) : Error opening SHP file

the working directory is in check and i also tried putting it inside the script to no avail.


texts: error running readLAS() in lidR: 'new_proj_and_gdal' is not an exported object from 'namespace:rgdal'

texts:
I've just updated lidR to 3.1.1 and am now unable to load a .las file using readLAS():
las <- lidR::readLAS(paste(sfm_dir,"/","201807XX_Denali_FHM_north_group1_densified_point_cloud_part_2.las",sep=""))
#> Error: 'new_proj_and_gdal' is not an exported object from 'namespace:rgdal'

I have not gotten this error with my previous version of lidR (3.0.1) which I just removed to allow for upgrade.
Maybe pertinent additional info:
R version 3.6.2 (2019-12-12)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)
rgdal_1.4-8
raster_3.3-13
sp_1.4-4


texts: I have a SpatialPointData Frame in R. Each point represents lat long and InsuredValue amount. I would like to find the maximum sum of insured value in a radius of 200 mt. How I can do in R


texts:
I'm trying to use Python to:

Do a union between 2 layers (LandTrust and Precincts) - this works
Create a field called "Overlap - this works
Create a da.UpdateCursor to populate the fields with values - this does Not work

I don't understand what I'm doing wrong (I'm using both the Python Window and an independent text editor).

Here is my code:

>>> import arcpy
... unionLayer = arcpy.Union_analysis(["LandTrust", "Precincts"], "UnionOutput", "ALL")
... arcpy.AddField_management(unionLayer, "Overlap", "TEXT", "", "", 20, "", "NULLABLE")
... updateRows = arcpy.da.UpdateCursor(unionLayer, ["FID_LandTrust", "FID_Precincts", "Overlap"])
... for row in updateRows:
...     if (row[0] != -1 and row[1] = -1):
...         row[2] = "LandTr"
...     elif (row[0] = -1 and row[1] != -1):
...         row[2] = "Pre"
...     elif (row[0] != -1 and row[1] != -1):
...         row[2] = "Overlap"
...     updateRows.updateRow(row)
... del updateRows
... del row

Parsing error SyntaxError: invalid syntax (line 6)
>>>

It seems that my if statement is the issue. I'm executing the code by pressing enter twice at the end.


texts:
I've have several orthophotos stored in an ECW (Enhanced Compression Wavelet) file format. For most of them the exact acquisition date is referred in the file name - e.g. "Paris_24_12_2011". But for some of them this information is lacking. I'm currently using ArcGIS 10.3.

Is there a way to check when the orthophoto was taken? maybe some metadata?
If so, how can I watch this data?

I've tried to look it up and found this page of ECW file specification including metadata for GDAL use. I'm not familiar with GDAL so I left it aside for the moment.
I've downloaded erdas er viewer 2014, thought it might give me some metadata access. It didn't.

Is there any simple, free-to-use way that can help me find the date of the orthophoto?


texts:
I have added delimited text layer from a csv with "watch file" option checked in and I want to update the mapcanvas every five second the program always crashes. I have attached the sample python script.

import time

while True:

    time.sleep(5)

    iface.mapCanvas().refreshAllLayers()


texts:
I'm trying to zoom to the extent of several vector layers data.

I've tried to create a new OpenLayers.Bounds bounds and to add a 'loadend' for each OpenLayers.Layer.Vector to increase bounds, but then I don't know where to put map.zoomToExtent(bounds)... is there a way to run this after all the layers have been loaded?

I've tried https://gis.stackexchange.com/a/6297, but it did not work, as the vector layers were not created yet...

I've also tried https://gis.stackexchange.com/a/41588/16356, but it did not work either, as I have several layers.

EDIT: Here is how I fixed the issue, thanks to Vadim.

bounds = new OpenLayers.Bounds();
map = new OpenLayers.Map('map');
LOADED_LAYERS = 0;

layers[1] = new OpenLayers.Layer.Vector("layer1", {
                eventListeners: {
                    'loadend': function (evt) {
                        LOADED_LAYERS++;
                        bounds.extend(layers[1].getDataExtent());
                        if (LOADED_LAYERS == 3) {
                            map.zoomToExtent(bounds);
                        };
                    }
                }
            });
layers[2] = ...;
layers[3] = ...;

map.addLayers(layers);
map.zoomIn();


texts: Exporting data frame as PNG in ModelBuilder using ArcMap

texts:
I've created a model that when finished creates two symbolized raster layers and adds them to the dataframe.
I would now like each layer to be exported as a unique png file.
Is there a tool that I can add to my model to complete this task?
See the model below. I would actually like to create two pngs, one of phaseFinished.lyr and another of activityFinished.lyr.


texts: Loading external WFS (GML output format) in Leaflet map

texts:
I need to build a Leaflet map loading data from this external WFS service

http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/wfs/Numeri_Civici_2012.map&service=WFS&request=GetCapabilities&

with these output formats

<OutputFormats>
   <Format>text/xml; subtype=gml/3.1.1</Format>
</OutputFormats>

Here you are a sample of the data

http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/wfs/Numeri_Civici_2012.map&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=IN.NUMERICIVICI.2012&SRSNAME=EPSG:4326&bbox=7.59,44.83,7.75,44.94&outputFormat=GML2&

I know how to load an external WFS if the output format: here you are a sample

  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.5/leaflet.css" />
  <title>Get External JSONP</title>

... but this obviously doesn't work with GML.

Any suggestions in case of GML output formats?


texts: How to assign node IDs to links in a network?

texts:
I have two shapefile layers in ArcGIS for Desktop 10.

I have one line layer and node layer.

My node layer(shp) has a "code_1" field and line layer has a "code_node" field.

The intersection point of the line layer has only one node.

I want get the code from the "code_1" field of the node layer for the line's "code_node" field

For example, for one line, the first node has code_1=1 and the node at the end of line has code_1=2. In this case the code_node attribute in that line must be "12",.

How can I automatically get the code from the nodes for each line?


texts: How to clear attribute columns from a Shapefile?

texts: I want to clear all columns from a SHP. Only need geometry of features. I use QGIS and Sextante with grass and saga. Is there any tool for this job...


texts: A QGIS tool to perform cost surface analysis (as in ArcMap Path Distance)? not a least cost path nor a network analysis?

texts:
Within QGIS I am trying to perform a surface analysis to show the part of the terrain that can be reached from one or more source points in let's say 1, 2 and 3 hours of walk on the ground. That is to determine the minimum accumulative travel cost from a source to each location on a raster surface (DEM), considering the actual surface distance as well as other horizontal and vertical factors (friction, slope, type of soil, etc.).
I am NOT looking for a least-cost path tool (that only considers two points), nor a network analysis tool (that uses a network of streets, and so on), and not even a simple buffer zone/Euclidean distance (that doesn't consider the slopes). I'd like to consider all the terrain around one or more sites, so that if the cost to travel between the points varied according to some characteristic of the area between them, then a given location might be closer, in terms of travel cost, to a different point.
I am attaching an image of a simple map I'd like to achieve that might help in understanding the issue.

I know that in ArcMap this task is done by the Path Distance tool (link), but I am not able to find any tool or plugin in QGIS, GRASS or SAGA. I found that the VARCOST in IDRISI should do the same (link).
The problem is also that I cannot find the exact term to define such analysis, as many different names seem to be used at the same time: cost surface analysis? path distance? travel time analysis? cost distance? travel cost? catchment area? service area? cost allocation? ...
I am using QGIS 3.16.1-Hannover on an Ubuntu 20.04.1 LTS laptop.


texts:
I have a problem that I don't understand. I was working with QGIS and Postgis on importing and exporting layer from and to database, everything was doing very well. Suddenly the connection to the database fails. I'm getting this error message in QGIS:

What could be the cause and how to solve it?


texts: Feature Services & arcpy.ConsolidateMap_management & .replaceDataSource

texts:
We have several in house published feature services, including some from our Esri Data Stores. When I use arcpy.ConsolidateMap_management on an mxd (10.1) the Feature Services are nicely consolidated into local file geodatabase, but the layers are not re-sourced.

When I use ArcCatalog to manually "set data source/s" in the mxd to the local file geodatabase, the paths are not present in the output mxd and no errors are produced.

Not stopping there. I then wrote some ArcPY to re-source the offending layers.The ArcPy works perfectly but the sources remain the same...

 if ISFS == 1: # Only proceed if it's an HTTP / Feature Service
            arcpy.AddMessage(" ")
            arcpy.AddMessage("     - FEATURE SERVICE TEST ---------------------" )
            arcpy.AddMessage("     - Layer Name: " + str(FSLyr))

        try:
            dataset_name = ("a"+str(FSFCNAME))
            arcpy.AddMessage("     - FS Data Source: " + str(LayerSource))
            arcpy.AddMessage("     - Consolidated fGDB: " + str(ConsolidatedFeatureServiceWorkspace))
            arcpy.AddMessage("     - Feature Class Name:" + dataset_name)
            FSLyr.replaceDataSource(ConsolidatedFeatureServiceWorkspace, "FILEGDB_WORKSPACE", dataset_name, False)
            arcpy.AddMessage("     - Replaced ")
            arcpy.AddMessage(" ")
            FeatureServiceMxd.save()

Results:

Before I investigate loading the file geodatabase feature class into the mxd and stealing the old cartography, does anyone know why we are denied changing Feature Service Layer sources?


texts: Where can I download sat imagery. I've read that OpenLayers provide such things but it requires an internet connexion, and I mostly use QGIS with an offline computer. So do you know a way/provider (for free) for such products?


texts: GPS drift during the receiver is stationed

texts:
I have developed a GPS Tracking device it is working fine when the receiver is in motion, when it stationed in a position for couple of hours i get some drift for around 10-15meters. The receiver send different lat/lon every 30seconds, which makes it impossible for me to calculate the total IDEAL time of the tracker device.

Is there any GPS accuracy filter which can be implemented in my software to overcome this issue?


texts: How to aggregate points to polygons using ArcGIS Desktop?

texts:
I am trying to aggregate my data points to the polygons so
that I could find out how many points I actually have in
each polygon, and then export the data contained within
each polygon to further look at the data properties as
they pertain to each polygon.


texts:
I am working with EPSG 27700 polygon data in PHP. All the polygons are in the UK. The polygons do not have holes inside them. I'm hoping to be able to calculate the area of these polygons. For this I believe I need to take two steps. Firstly to convert the coords into something where x an y are the same distance and then once that is done to calculate the area inside each polygon.

My naive solution would be to use some graphics package to draw the polygon onto an image and count the coloured squares.

What's a good coordinate system for this? UK Easting & Northing is my best guess.
Is there a PHP library or easily(!) implemented algorithm to find the area inside a polygon defined by a list of points?

Alternatively, should I be getting out of my comfort zone and learning to use a proper GIS system?


texts:
I'm writing a web app using ESRI's Javascript API and a map I created in ArcMap 10.1 hosted on ArcGIS online (Feature Service, not WebMap).

I have written the code necessary to print attributes in a stationary DIV instead of creating a popup. The only problem is, I now have to click twice every time I want to select a new feature.

So, I select a feature, the DIV updates with the new attributes. I click on a different feature, nothing happens. I click on the same feature again and the attributes update in the DIV.

What is going on here? I'm using map.infoWindow.on() to listen for "selection-change". My code is below and the website is live here.

popup = new Popup({popupWindow:false},domConstruct.create("div"));

map = new Map("Div1", {
 basemap: "gray",
 center: [-86.14, 39.79],
 zoom: 11,
 infoWindow: popup,
});

featureLayer = new FeatureLayer("http://services1.arcgis.com/ROSiDsWkQXhsUAs0/arcgis/rest/services/Thesis.gdb-2015-01-19_(1)/FeatureServer/0",{
 mode: FeatureLayer.MODE_ONDEMAND,
 infoTemplate: popupTemplate,
 outFields: ["*"]
});

map.addLayer(featureLayer);

popup.on("selection-change", function(){
 var feature = popup.getSelectedFeature();
 var content_MED_dolla = feature.attributes["MED_dolla"];
 document.getElementById("popupDiv").innerHTML = "
Median mortgage value: $" + content_MED_dolla + ",000
";

});


texts: PyQGIS/QGIS2.18 : Store checkbox state as boolean in database (PostGreSQL 9.3)

texts:
I would like to store the information about whether a checkbox is checked or not in a boolean column in my database (running on PostGreSQL 9.3). I'm using a custom form, and the corresponding checkbox field should appear as checked or unchecked depending on its value in the database - and the value in the database should update if the user changes the state of the checkbox.

I created a variable for the checkbox in my Python script:

global zone_humide
zone_humide = dialog.findChild(QCheckBox, "zone_humide")

I also made sure that the field has the same name as its corresponding column in the database.
But when I check or uncheck the box, the change is not recorded in the database - and consequently, I can't call back the current state of the box when I open my form again.