texts: Reproject MODIS data using R (results in NAs or no spatial extent)

texts:
I am using GLASS albedo data stored here for pre-2000 (AVHRR) data and here for post-2000 data (MODIS). My end goal is to create a raster stack of each month that contains white sky albedo data from 1982-2015. The problem I have run into is that the MODIS and AVHRR data are in different spatial reference systems and I can't seem to reproject them to be in the same system.

I convert from hdf to tif using R like this:

fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
        ".../modis.tif")
gdal_translate(get_subdatasets(fileavhrr)[8], projwin = c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like data north of 50 degrees

avhrr<- raster(".../avhrr.tif")

#class       : RasterLayer
#dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
#resolution  : 0.05, 0.05  (x, y)
#extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

modis<- raster(".../modis.tif")

#class       : RasterLayer
#dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
#resolution  : 154.4376, 308.8751  (x, y)
#extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
    +b=6371007.181 +units=m +no_defs
#values      : -32768, 32767  (min, max)

Here are things I have tried:

1.) Use the MODIS Reprojection Tool. For whatever reason, this tool seems to think the subdatasets of the MODIS .hdf files are only one tile (the upper left most tile, tile 0,0) and not the global dataset. My understanding is that the MODIS data are global (not in tiles?), so I do not know why the MRT is doing this.

2.) Use the raster package in R.

projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")

This returns a raster with values that are all NA:

class       : RasterLayer
dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
values      : NA, NA  (min, max)

3.) Use the gdalUtils package in R:

gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile= ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )

This returns a raster with essentially no spatial extent.

gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
#class       : RasterLayer
#dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
#resolution  : 0.02801551, 0.02801573  (x, y)
#extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

Any ideas on why reprojecting this MODIS data is so difficult?


texts: Converting a table to a raster for clipping

texts:
I have a series of tables (TAB files) that represent the calculated values for a corresponding .tif image.  These tab files have the same number of rows and columns as pixels in the tif images.  The tif files are georeferenced and correspond to a study area I am using and I want to be able to extract the tab info for only the values within the study area.  Anyone have any suggestions on how to go about this?

Or if I can just relate the table to the tif so each column and row is related to the corresponding pixel of each column/row in the tif.

I'm using ArcGIS 10.1 for Desktop.


texts:
I've several files composed of elements like the one below which I'm trying to make sense of (essentially a very simple microtopography study).

I've vectorised the image with lots of separate lines ending where they intersect with another line. Each of these vertices has been assigned a separate number corresponding to a line, and these lines are in turn arranged into units (so the photo below is of 1 unit composed of 8+ lines composed of multiple vertices).

I'd like to analyse these to look for regularity in the number of lines comprising a unit, does anyone know of doing this? I would use the "Count unique features in polygon" only I'm only using lines.

I'm using QGIS 2.18 Las Palmas


texts: Finding all the POIs in the world using OpenStreetMap?

texts:
I would like to extract airports' boundaries from OpenStreetMap.
Therefore, to extract data from a given country, for instance France I do:

area["name"="France"]->.boundaryarea;
(
node(area.boundaryarea)["aeroway"="aerodrome"];
way(area.boundaryarea)["aeroway"="aerodrome"];>;
relation(area.boundaryarea)["aeroway"="aerodrome"];>>;
);
out meta;

However, I would like to search for such data for the whole world, i.e. without limiting the area.
Is there a node representing the world or do I have to loop such queries country by country ?


texts: Searching with multiple parameters in ArcGIS REST API?

texts:
I need to perform a search on a MapServer that populates a datagrid with the resulting information such as in this example Show find task results in a DataGrid. However, I need to be able to have multiple parameters to search for at the same time. For instance, I need to search for a person at a specific address.

I am using findParameters() and findTask() but I don't see how I can use multiple parameters in a searchText.


texts: Clipping Raster Mosiac Simultaneously Within Separate ArcPy Scripts?

texts:
I have a large raster mosaic dataset that I would like to clip within several scripts simultaneously, but whenever I run two instances of a script simultaneously one always fails with the error:

  arcgisscripting.ExecuteError: ERROR 999999: Something unexpected
  caused the tool to fail. Please refer to the error help for potential
  solutions, or contact Esri Technical Support
  http://esriurl.com/support. Failed to delete raster dataset Failed to
  execute (Clip)

Is there a way to allow multiple scripts to clip a mosaic raster dataset at the same time? Here is some example code:

import arcpy
import os
ras = arcpy.Raster(os.path.join('d.dbf', 'mosaic'))
pol = 'pol.shp'
result = arcpy.Clip_management(ras, '', 'result.tif', pol, '', 'ClippingGeometry')


texts: ArcPy Insert cursor losing vertices when writing Polygon features?

texts:
I'm currently trying to code a geodesic buffer from points (the one that is implemented in ArcGIS is very slow). The idea is to calculate 100 point coordinates (P_Buffer) around an initial point (P) and to shape a polygon from the P_Buffer Points (accordingly, the buffer polygon has 100 vertices).

The script is structured like this:

Create Polygon Shapefile
Open input file and get coordinates for P
Loop: Calculate P_Buffer coordinates (Vincenty's formulae) and create an arcpy Point Geometry (100 times)
Create an arcpy Polygon Geometry from 100 Point Geometries
Write Polygon Geometry (Insert Cursor)

That's the script (simplified):

# Step 1

infile = "CRATER_Pickering.shp"
outfile = "C:\test.shp"
sr = arcpy.Describe(infile).spatialReference
arcpy.CreateFeatureclass_management("C:/", "test.shp", "POLYGON", "", "", "", sr)
cur = arcpy.da.InsertCursor(outfile, ["SHAPE@"])

# Step 2 & 3 work perfectly fine - I skip most part of the section. Each time the loop runs,
# step 3 ends with adding the P_Buffer coordinates to an arcpy Point Array. After that,
# the arcpy Point array is appended to an arcpy Array, which stores all the P_Buffer
# coordinates / the buffer polygon's vertices.

...(get buffer distance and coordinates of P)...

    polygon = arcpy.Array()

        ...(coordinate transformation)...

        point = arcpy.Point(P_Buffer_X, P_Buffer_Y)
        polygon.append(point)

# Step 4

    polygon_shape = arcpy.Polygon(polygon)

# Step 5

    cur.insertRow([polygon_shape])

Everything works fine until I get to step 5. The coordinates are calculated correctly (just as in the ArcGIS geodesic buffer) and all 100 P_Buffer points are polygon_shape vertices in Step 4. However, when writing the geometry in step 5, only few buffer polygons have 100 vertices. Big circular buffer polygons have 100 vertices but the smaller their diameter gets, the fewer vertices the polygons have (up to zero - meaning that smaller buffer polygons either show angular boundaries or that no geometry is visible).

Here's an example. The green dots in the image show the P_Buffer points from the polygon array in step 3. The buffer polygons written in step 5 are displayed in pinkish colour.

Huge diameters are fine:

Small diameters are not:

How is this possible? When I write a Point Shapefile from the polygon array  (Step 3 - all 100 vertices are still there) or count the vertices in polygon_shape (step 4 - all 100 vertices are still there), everything looks fine. As soon as I call insertRow in Step 5, vertices are excluded. Do you have any idea what I did wrong?

Solved:

It was indeed a spatial reference issue. As suggested by klewis, I added the spatial reference to step 4. Now the polygons have all the vertices from the coordinate transformation included.


texts: Missing date column from interactivity callback

texts:
I can't seem to get "date" columns to return for popups/infoWindows. I can retrieve other data types without an issue. I have the interactivity set correctly from what I can tell, but the featureClick callback only returns one column.

My code relies on v2 of the CartoDB API.

Here's my call to cartoDB.createLayer:

cartodb.createLayer( map, { type : "cartodb", options : layerData } )

where layerData is

interactivity: "name,admin_date,size"
query: "SELECT size, name, admin_date, the_geom_webmercator FROM book1"
table_name: "book1"
tile_style: "#book1{marker-width:5;marker-opacity: 0.9;marker-line-color: white;marker-fill:#e6e6e6;marker-line-width: 1;marker-line-opacity: 0.9;marker-placement: point;marker-type: ellipse;marker-allow-overlap: true;} #book1{[size>=0]{marker-fill:#b8cde9;}[size>=18229430]{marker-fill:#6693cc;}[size>=22061458]{marker-fill:#2f5f9d;}}"
user_name: "c2es"

The featureClick callback only return the first and last columns. Any ideas what's wrong?


texts: Exporting list of lists in Google Earth Engine

texts:
I have a list of lists (you can see on the pic), and would like to export it as .csv. I tried to convert it into feature collection, and then run the following lines:
var feature = ee.Feature(null, lossperyearperelevation1);
var featureCollection = ee.FeatureCollection([feature]);
Export.table.toDrive({
  collection: featureCollection,
  description: 'exportTableExample',
  fileFormat: 'CSV'});

it gives the following error: Error: Feature, argument 'metadata': Invalid type. Expected type: Dictionary. Actual type: List<List<Dictionary>>.
I am confused if it is a list of list, and how to export it in .csv. How can I solve this problem?
the link for the sample code is here: https://code.earthengine.google.com/40c2c24233a6119ee577fb05dd86a14e


texts: Export Geotif with transparent Background

texts:
i using Arcmap and i have a raster file, i need export this using a Geotiff file for use in Geoserver like a layer.

My problem is, when i export the file always exist a color gray or white in the background.

In frame properties i select white background with 'no color', when export geotiff in format tab, i select in background color 'no color'.

Exist another way for found some solution?


texts: Connecting table and feature class in ArcGIS for Desktop?

texts:
I have a polygon feature class inside of my geodatabase, and I want to attach attributes from some tables that resides in a separate database. Is there any way to do this that would maintain the relationship?

I have it working with joins, but that doesn't keep the relationship inside of the feature class.


texts: How to show shapefiles with Bing map in asp.net

texts:
I want to display my shapefile data on bing map using ASP.NET.

Its an web application.

I don't want to use any third party tool for shapefile conversion.


texts: Use high precision to convert image from one EPSG defined CRS to another

texts:
This question is based on another question: I'm converting multiple images (.jpg or .png) from one EPSG defined CRS system to another, e.g. from EPSG:31255 to EPSG:25833 (on a server within a Java app).

Currently I'm using this code (GeoTools 22.2, Java 8):

//codeIn file (already includes TOWGS84): http://epsg.io/31255.wkt
String codeOut = "epsg:25833";

public void convert(String pathIn, String pathOut, String codeOut) {
    File fileIn = new File(pathIn);

    if(fileIn.exists() && fileIn.isFile()) {
        try {
            File fileOut = new File(pathOut);
            CoordinateReferenceSystem crsOut = CRS.decode(codeOut);
            AbstractGridFormat format = GridFormatFinder.findFormat(fileIn);
            Hints hints = null;

            AbstractGridCoverage2DReader reader = format.getReader(fileIn, hints);
            GridCoverage2D coverage = reader.read(null);
            reader.dispose();

            GridCoverage2D coverageTransf = (GridCoverage2D) Operations.DEFAULT.resample(coverage, crsOut);
            GeoTiffFormat outFormat = new GeoTiffFormat();
            GridCoverageWriter writer = outFormat.getWriter(fileOut, hints);
            writer.write(coverageTransf, null);
            writer.dispose();
        } catch(Exception e) {
            e.printStackTrace();
        }
    }
}

Now I want to add high precision to it. In Proj4 I'd use this command:

+ellps=bessel +nadgrids=GIS_GRID.gsb +proj=tmerc +lat_0=0 +lon_0=13d20 +k=1 +x_0=0 +y_0=-5000000 +units=m +no_defs +geoidgrids=foo.gtx +to +datum=WGS84 +proj=utm +zone=33 +no_defs

In this similar question there's a link to a website that tells you to copy the two files (.gsb and .gtx) to the user_projections folder, which of course doesn't exist in my app.

My guess is that I'll have to give GridCoverage2D coverage = reader.read(null) some type of GeneralParameterValue[]. Is this the right way to do it and if so, how exactly do I give it the path of the two files and/or the command?


texts: Polygon Perimeter from Raster with Conditions

texts:
How can a polygon be created around a group of raster cells in ArcMap?

Specifically, given the set of red cells in the figure, I'd like to create a polygon that represents the perimeter. The criteria is that the perimeter excludes cells which are disconnected from the main red area by one or two white cells. I do want to avoid coding in ArcMap.


texts: Automatically updating Buffer FeatureClass when creating new features

texts:
I have a roads FeatureClass, and I've built a buffer for it. I will still make edits to the road FeatureClass.

Let's say I add a new segment of road, how can I automatically buffer the new segment without creating a new buffer FeatureClass?

I use ArcGIS Desktop 10.1

Apparently when using Geomedia if your information changes in your data source, so does the analysis in the workspace. I want to know if ArcMap can do the same.


texts: Conflicting measurement units with buffer ftools

texts:
I am new to qgis and I would like to create a buffer around point features.

My goal is to create a quarter mile buffer (402.336m) around each point feature, using Geoprocessing Tools -> Buffer(s), setting the buffer distance and exporting. The process creates what appears to be a massive polygon. However when I use '.05' as the distance buffers appear around each point.

I assumed this happened because the measurement units are not correct. I verified the project projection (NAD83 OH/South). Then I checked the Measure Tool via Options -> Map Tools. Te ellipsoid for distance calculations set at WGS 84 and the Preferred Measurement Units is set to 'meters'

I have set up 'on the fly' projections

How do I solve the mystery of this buffer problem?

Setup: qgis 1.7, Ubuntu 11.04


texts: Select by Attribute --Top 500 values selection in attribute table - ArcGis

texts:
Is there any way we can select top 500 values in attribute table?
I tried by sorting area field ascending/descending... but I can't count big numbers to select.
In microsoft excell we can do know, similarly in arcgis, is there any code to keep in attribute selection ??


texts: Turn Layer off by default in Leaflet LayerControl

texts:
In the following python code I have three data layers while I need to show only two of them by default as the program executes. So the Places and District layer should be checked as the map is opened and the Country layer should be unchecked/turned off. What should I write in LayerControl for this?

import folium
import pandas

data=pandas.read_csv("Previous.csv")
data = data[data['Death'] > 0]
lat = list(data["Latitude"])
lon = list(data["Longitude"])
dat = list(data["Date"])
are = list(data["Area"])
sym = list(data["Symptoms"])
tra = list(data["Area_traveled"])
sou = list(data["Source"])
dea = list(data["Death"])

def color_producer(elevation):
    if elevation <= 1:
        return 'orange'
    elif 1 < elevation <= 3:
        return 'red'
    elif 3 < elevation :
        return 'darkred'

fga = folium.FeatureGroup(name="Places")
map = folium.Map(location=[23.789412, 90.351756], zoom_start=7, tiles="Mapbox Bright")
for lt, ln, co, mo, vi, up, na, de  in zip(lat, lon, dat, are, sym, tra, sou, dea):
    fga.add_child(folium.CircleMarker(location=[lt, ln], radius = str(de), popup= "
Link " + na + "
Area: " + co + "
Symptoms : " + vi + "
Area traveled : " + mo+ "
Source : " + up,
    fill = 'true', fill_color = color_producer(de), color = 'grey', fill_opacity=0.7))

fgp = folium.FeatureGroup(name="District")
fgp.add_child(folium.GeoJson(open("dist.geojson",encoding = "utf-8-sig").read(), style_function=lambda x: {'fillColor':'grey' if x['properties']['Death'] == 0
else 'darkred'}))

fgq = folium.FeatureGroup(name="Country")

fgq.add_child(folium.GeoJson(open("world.json",encoding = "utf-8-sig").read(), style_function=lambda x: {
'fillColor': 'none'}))

map.add_child(fga)
map.add_child(fgp)
map.add_child(fgq)
map.add_child(folium.LayerControl())
map.save("Corona.html")


texts: How to efficiently allocate memory for big raster file in lidR v2.2.4?

texts:
I am intended to find trees from a CHM created in agisoft with total raster size of 1.2GB with this following code:

library(lidR)
library(future)

# set to use 8 cores with future
plan(multisession, workers=8L)

mainDir <- "E:/02.CHM/01.Agisoft_Analysis/South"
inDir <- "4_ground_classified"
outDir1 <- "5_trees"
rasterfile <- "CHM_0.1_South.tif"
treefile = "tree_south.csv"

dir.create(file.path(mainDir, outDir1), showWarnings = TRUE)

chm <- raster(file.path(mainDir, inDir, rasterfile))
ttops <- tree_detection(chm, lmf(ws=5, hmin=1.37, shape=c("circular") ))

When I was processing the raster with tree_detection() function in lidR v2.2.4, it failed with this message:

Error: cannot allocate vector of size 647.2 Mb.

I tweaked by creating a tiled raster with this command

chm <- raster(file.path(mainDir, outDir1, "CHM_0.1_South.vrt")
ttops <- tree_detection(chm, lmf(ws=5, hmin=1.37, shape=c("circular") ))

But the same error message still appears.

The second tweak was by using memory.allocation() trick, but that also didn't work.

As for your information, my RAM size is 16GB.

I found this vignette to work on parallel computation for las files, but I can't find a way to do parallel computation on VRT .tif file.

Do you have any suggestion on how to do parallel processing computation for tiled .tif raster?
Given the fact that R can't allocate enough memory if I work in one big .tif raster file.


texts: Why do I get an error running Catchment Polygon Processing?

texts:
I am trying to run catchment polygon processing in ArcGIS 10.3, and this is what i get "error hresult e_fail has been returned from a call to a com component"
Can you help me to debug this error?


texts: How to extract data from planet_osm_ways and insert into planet_osm_polygon again?

texts:
I have an OSM import in my PostgreSQL database and can see my map using Geoserver.

I've noticed some data have disapeared ( most water related data ) maybe after an import failure, I don't know.

I'm thinking in delete all water from database and then import again from the ways table because I believe all data still in there to allow append imports.

Another alternative is to drop the database and do a full import again but it will be a pain because it is a planet database (almost two weeks to do the import).

Some ideas ?


texts:
I'm trying to add the pixel coordinates (latitude and longitude) as bands to export (code below). However, the export returns 1 decimal place for the latitude and longitude. How can I increase the precision (# decimals) in my export for the coordinates?

//function to add coordinates
var getLatLong = (function(img) {

  var proj = img.select([0]).projection()
  // get coordinates image
 var latlon = ee.Image.pixelLonLat().reproject(proj)
 var img1 = img.addBands(latlon.select('longitude','latitude'))

  return img1;
});


texts: Loading GML data into Postgres - No geom!

texts:
I'm trying to load an Environment Agency GML file into Postgres/PostGIS using ogr2ogr. The file loads fine, all the other data is there, but when I check the table in the database there is no geometry column. As it is a polygon dataset this is a problem!

ogr2ogr command is:

ogr2ogr -f "PostgreSQL" PG:"host=* port=* dbname=* user=* password=*" RoFRS_v201701.gml -overwrite -progress --config PG_USE_COPY YES

The GML file has a XSD file with it.

Also the GML opens perfectly in QGIS, with polygons displayed as they should be.

Any suggestions for getting the geom into Postgres?

Update:

Results from ogrinfo:

Geometry: Polygon
Feature Count: 2130752
Extent: (134111.801000, 11510.471300) - (655988.124900, 657600.687500)
Layer SRS WKT:
(unknown)
OBJECTID: Integer (0.0)
PROB_4BAND: String (8.0)
SUITABILITY: String (25.0)
PUB_DATE: Integer (0.0)
RISK_FOR_INSURANCE_SOP: String (3.0)
SHAPE_Length: Real (0.0)
SHAPE_Area: Real (0.0)


texts: Show only one label for multiple points with same value in QGIS

texts:
I have a point shapefile which has for a few occasions multiple points with the same ID. Doesn't look nice with the normal label-display.
Is it possible in QGIS to only display one label then? I couldn't find a check/uncheck option... Perhaps it would be possible to retrieve labels via a SQL-query, with a »GROUP BY id«?


texts: Can Geoserver let users choose a single tile from a layer to be displayed?

texts: I have a raster layer of about 30 tiles. The rasters overlap and important data is hidden when the tiles are displayed as a layer. I'm wondered how others get around this issue? Can geoserver let users switch on and off individual tiles in a raster layer?


texts: Disable doubleclick zooming in OpenLayers 3

texts:
I've searched for it on the entire web and finally I found on artamstrong.com how to disable a single feature like "mouseWheelZoom".

Here is how doubleclick zooming in OpenLayers 3 can be disabled:

var map = new ol.Map({
    controls : ol.control.defaults()
                .extend([ new ol.control.FullScreen() ]),
    interactions : ol.interaction.defaults({doubleClickZoom :false}),
    target : 'map',
    layers : [ new ol.layer.Tile({
            title : 'OpenStreetMaps',
            preload : Infinity,
            source : new ol.source.OSM(),
            visible : true
    }) ],
    view : new ol.View({
            center : ol.proj.transform([ 9.41, 48.82 ], 'EPSG:4326','EPSG:3857'),
            zoom : 12
    })
});

It's really that simple, just tell ol3 to set "doubleClickZoom" on false:

interactions : ol.interaction.defaults({doubleClickZoom :false})


texts: How can I tell what licensing level my ESRI GIS server is from client tools?

texts:
From client tools, like ArcCatalog, is there a way to know if the GIS server I am connecting to is running under a workgroup or enterprise license?  Also, is it possible to tell if the edition is basic, standard, or advanced?

I ask because we are trying to migrate a 9.3.1 database to 10 (migrate but not upgrade yet) using ArcCatalog.  When we add a new database server in ArcCatalog the only options are to add a personal or workgroup database server.  It seems like it is only looking for SQL Express instances, but it should be also be able to see full SQL Server (Enterprise) servers.  We have ESRI enterprise licensing, but I want to know if it is possible to discover that from the client tools.

Help, About ArcCatalog shows ESRI ArcCatalog 10, ArcGIS Desktop 10 SP2, license type ArcInfo.  This is the ArcCatalog (client) licensing, not the server.


texts: GDAL Python Function to convert vector files from shapefile to KML

texts:
I want to convert a shapefile with polygon to KML file using GDAL Python.
For this I created a fuction with two arguments (input and output) and it works, but points in polygon are switched (X in KML is Y in SHP and Y in KML is X in SHP).
In ogr.ogeometry is a method SwapXY but it is for points not polygons.
https://gdal.org/python/osgeo.ogr.Geometry-class.html#SwapXY
I try to export polygon to array of points but I didn't find any fuction for it.
After that I create geometry from WKT and export to KML but I don't know how to save it into output file (not as text file but KML with georeferences).
I don't want use subprocess with ogr2ogr, because it is a system dependency.
Here is a code:
    @staticmethod
    def __convert_shp_to_kml(shp_file: pathlib.Path, kml_file: pathlib.Path = None):
        """
        Convert shp (shape file) to kml file.

        :param shp_file: path to shp file
        :type shp_file: pathlib.Path
        :return: True if conversion is correct
        :rtype: bool
        """
        assert shp_file.exists() is True, logging.error("File: " + str(shp_file) + " not exist.")
        if kml_file is None:
            kml_file = str(shp_file).replace(shp_file.suffix, ".kml")
            kml_file_name = str(shp_file.stem)
        else:
            kml_file_name = str(kml_file.stem)
            kml_file = str(kml_file)

        driver = ogr.GetDriverByName('ESRI Shapefile')
        dataset = driver.Open(str(shp_file))
        layer = dataset.GetLayer()
        spatial_ref = layer.GetSpatialRef()
        epsg_code: int = int(spatial_ref.GetAttrValue('AUTHORITY', 1))

        output_driver = ogr.GetDriverByName('KML')
        out_ds = output_driver.CreateDataSource(kml_file)
        output_data_set = osr.SpatialReference()
        output_data_set.ImportFromEPSG(epsg_code)
        out_layer = out_ds.CreateLayer(kml_file_name, output_data_set, geom_type=ogr.wkbPolygon)

        ds = ogr.Open(str(shp_file))
        if ds is None:
            logging.warning(str(shp_file) + "ds is empty")
        lyr = ds.GetLayer()
        for feat in lyr:
            geom_poly = ogr.CreateGeometryFromWkt(str(feat.GetGeometryRef()))
            print(geom_poly.ExportToKML())
            out_feat = ogr.Feature(out_layer.GetLayerDefn())
            out_feat.SetGeometry(geom_poly)
            out_layer.CreateFeature(out_feat)
            out_layer.SyncToDisk()
        return True


texts:
I'm trying to add vertical and horizontal rulers on top of my map. In OpenLayers 2 I have seen code that had VerticalRuler and HorizontalRuler as attributes but I can't find something similar in version 6. Has it been removed as a feature or has it been renamed?

An image that showcases how it should look:


texts: Speeding up connection to MS SQL Server in QGIS?

texts:
Currently when I try to connect to a SQL database in QGIS it takes a bout 30 seconds for all features in the database to load.

I have over 100 tables with geometry columns. When I try to connect to sql server it takes some time to read the table geometries and srid from SQL Server. I have seen the implementation of a geometry_column table in SQL Server for QGIS before. Is this still the most effective approach?


texts: What licenses are needed for ArcGIS Server and SDE (vs a FOSS setup)?

texts:
Does anyone know where I could find a cost breakdown for an ESRI/Microsoft GIS server and database setup?

I'm trying to figure out all the costs associated with an ArcGIS Server/SDE sitting on top of SQL Server, and compare it to the support cost of a FOSS setup.

From Microsoft we appear to need:

Windows Server - cost per processor and CAL for each user (ArcGIS Server apparently needs IIS which is part of Windows Server)

SQL Server - cost per core and CAL for each user

From ESRI:

ArcGIS for Server, Server Image Extension, ArcGIS for Desktop Standard Level (Needed to edit data in SDE)

Any information is greatly appreciated.


texts: Can be Neo4j Spatial data exposed to GeoServer via the Neo4j Geoserver Plugin being edited via WFS-T protocol?

texts:
I'd like to set up some data storage as Neo4j as I need the advantages of graph databases for may data, but I'm not sure that the data would be spatially editable via WFS-T and Geoserver. Does Neo4j Geoserver plugin just allows showing the data, or it make possible the spatial edit too?

Is that is not the case, that is the technologies stack you could recommend to enable the spatial data change for the Neo4j spatial database?


texts:
I am using three cons with raster calculator to give four different results depending on the values within two rasters, aclip20151111 and sclip20151111. There are four results I want; RESULTA, RESULTB, RESULTC AND RESULTZ, RESULTZ being the result for when none of the three conditions are met. My code is below:

CON(((0 <= "aclip20151111" <= 45 | “aclip20151111” > 315)&(“sclip20151111” <
3)),RESULTA,CON(((45 < "aclip20151111" <= 135 | 225 < “aclip20151111” >= 315)
&(“sclip20151111” < 10)),RESULTB,CON(((135 <= "aclip20151111" <= 225)&(“sclip20151111” < 15)
),RESULTC,RESULTZ)))

I am getting the error "ERROR 000989
Python syntax error: Parsing error SyntaxError: invalid syntax (line 1)".

What is the error in my input to raster calculator? I am using raster calculator as a tool in Arcmap.


texts: Creating new Feature Layer from existing feature layer on ArcGIS Online using ArcGIS API for Python

texts:
I have a folder on ArcGIS Online for a city that has few layers in it.
I need to create separate folders for each city and then create new feature layers using these existing layers, preferably using Python script.
I'm stuck right as I don't know how to create new feature layer from existing layer using Python (also, this new layer shouldn't be linked to the existing layer. The configuration should be exact same with only difference that they should have separate id's).
from arcgis.gis import GIS

gis = GIS("https://www.arcgis.com", username, password)
# Log into ArcGIS Online
gis = GIS('home')

#Town input from user
town_name = input('Enter town name: ')
town_abr = input('Enter town initials/Site id: ')

#creating town folder
gis.content.create_folder(folder='STR_{}'.format(town_name))

search_result = gis.content.search(query="title:STR_ENK1_COL_Duct Route Existing", item_type="Feature Layer")
(search_result[0]).id

existing_layer = gis.content.get('f21d1f7edbce478e88782e5a29728c11')


texts: Converting coordinate and raster values of each pixel to csv of large raster using Python?

texts:
There should be a simple way to convert a large raster to a csv or dataframe by reading the raster in as chunks, converting raster chunk to table of just three columns (x, y, and values), and outputting to csv or dataframe. I normally use something like the following in R but it only works for rasters much smaller than I'd like and I'd prefer to use python.

ras = readGDAL("Agroforestry_crop_setnull.tif")
rasdf= as.data.frame(ras, row.names=NULL, optional=FALSE, xy=FALSE, na.rm=TRUE)

I have read about GDALs ReadAsArray and gdal2xyz, I do not know how to use ReadAsArray to get coordinates of pixels in the 3 column table I desire, although it does allow reading the rasters as chunks. gdal2xyz does not provide the raster values of the pixels and does not allow reading as chunks from what I have seen.

In the past I have made use of ArcPy, breaking up the rasters into pieces, converting to points, and adding the coordinates to the attribute fields, then exporting these tables. I am hoping there is a better way, similar to the R solution I posted but for python and big rasters.

Is there a suggestion on reading big rasters as chunks in the csv format I desire -3 columns (X,Y, value)?


texts: Points in Polyon with two different CRS using QGIS?

texts: I'm relatively new to QGIS. I've got a layer that shows the boundaries of the wards in the UK. In addition to this, I've got a csv file, which has long/lat plus an answer (binary). The problem is, that these two files use different CRS (the one with the wards uses a "user defined coordinate system (USER 100001), and the point layer WGS 84. If I try to change them, either the points disappear or the wards... I actually wanted to find out how many points there are in one Query, but this doesn't work.. I assume because of the different CRS.. :( Is there a way I can change this? I've attached a screenshot.


texts: Error creating a spatial database. ERROR : could not load library "/usr/pgsql-9.1/lib/rtpostgis-2.0.so"

texts:
I am running PostgreSQL 9.1 with Postgis 2.0 on Fedora 15.

When trying to install the raster support,

psql -d mydb -f rtpostgis.sql -v ON_ERROR_STOP=1

I get the following error

psql:rtpostgis.sql:46: ERROR:
     could not load library "/usr/pgsql-9.1/lib/rtpostgis-2.0.so": libhdf5.so.6:
     cannot open shared object file: No such file or directory

A locate libhdf5.so.6 gives the following paths.

/usr/lib64/mpich2/lib/libhdf5.so.6
/usr/lib64/mpich2/lib/libhdf5.so.6.0.4

Reinstalling the GDAL from the repositories did not solve the problem. The version of the installed GDAL.

gdal.x86_64                                1.7.3-14.fc15                @updates


texts: Vender 0ption autoWrap is not working Geo Server 2.10.1

texts:
I have added a vender option of code for autoWrap to my xml and it does not seem to be working.  I have used multiple different pixel lengths to see if I could see the update but I have not been able to see any change as a yet.  In my reading, I came across a statement that autoWraping needs to be enabled.  as I went further down the rabbit hole I could not find any documentation on that.  listed below is the screenshot of the polygon I am watching to see if the .xml changes take as well as a snippet of the code I am editing.

         <sld:Rule>
        <ogc:Filter>
          <ogc:PropertyIsLike wildCard="*" singleChar="#" escapeChar="!">
            <ogc:PropertyName>building</ogc:PropertyName>
            <ogc:Literal>*</ogc:Literal>
          </ogc:PropertyIsLike>
        </ogc:Filter>
        <sld:MinScaleDenominator>0</sld:MinScaleDenominator>
        <sld:MaxScaleDenominator>3600</sld:MaxScaleDenominator>
        <sld:Name>Yes</sld:Name>
        <sld:Title>Building</sld:Title>
        <sld:PolygonSymbolizer>
          <sld:Fill>
            <sld:CssParameter name="fill">#d4d6e4</sld:CssParameter>
          </sld:Fill>
          <sld:Stroke>
            <sld:CssParameter name="stroke">#d4d6e4</sld:CssParameter>
            <sld:CssParameter name="stroke-width">0.5</sld:CssParameter>
            <sld:CssParameter name="stroke-linejoin">round</sld:CssParameter>
          </sld:Stroke>
        </sld:PolygonSymbolizer>
        <sld:TextSymbolizer>
          <sld:Label>
            <ogc:PropertyName>name</ogc:PropertyName>
          </sld:Label>
          <sld:Font>
            <sld:CssParameter name="font-size">12</sld:CssParameter>
            <sld:CssParameter name="font-family">Nimbus Sans L Regular Condensed</sld:CssParameter>
          </sld:Font>
          <sld:LabelPlacement>
            <sld:LinePlacement>
              <sld:PerpendicularOffset>0</sld:PerpendicularOffset>
            </sld:LinePlacement>
          </sld:LabelPlacement>
          <sld:Fill>
            <sld:CssParameter name="fill">#749dae</sld:CssParameter>
          </sld:Fill>
          <sld:VenderOption name="autoWrap">20</sld:VendorOption>
          <sld:VendorOption name="maxDisplacement">50</sld:VendorOption>
          <sld:VendorOption name="labelAllGroup">TRUE</sld:VendorOption>
          <sld:VendorOption name="removeOverlaps">TRUE</sld:VendorOption>
          <sld:VendorOption name="followLine">TRUE</sld:VendorOption>
          <sld:VendorOption name="group">TRUE</sld:VendorOption>
          <sld:Halo>
            <sld:Radius>
              <ogc:Literal>0.26</ogc:Literal>
            </sld:Radius>
            <sld:Fill>
              <sld:CssParameter name="fill">#bce6f4</sld:CssParameter>
            </sld:Fill>
          </sld:Halo>
        </sld:TextSymbolizer>
      </sld:Rule>


texts: Remove many menubar in startup.py on QGIS

texts:
I want to remove (not all) menubar with python. I will put the code in startup.py so that when loading QGIS, automatically just only one menu that is 'Project' menu. I get the code from Finding the children name of menubar in QGIS Python.

@xunilk and @Joseph gave the answer. The code :

from PyQt4.QtGui import QToolBar, QDockWidget, QMenuBar
from qgis.utils import iface

for x in iface.mainWindow().findChildren(QMenuBar):
    print x.objectName()
for item in x.actions():
    print item
for i, item in enumerate(x.actions()):
    print i+1, 'name: ', item.icon().name()
for item in x.actions():
    if item.text().replace('&','') == 'Edit':
       item.setVisible(False)
    if item.text().replace('&','') == 'View':
       item.setVisible(False)
    if item.text().replace('&','') == 'Layer':
       item.setVisible(False)
    if item.text().replace('&','') == 'Settings':
       item.setVisible(False)
    if item.text().replace('&','') == 'Plugins':
       item.setVisible(False)
    if item.text().replace('&','') == 'Vector':
       item.setVisible(False)
    if item.text().replace('&','') == 'Raster':
       item.setVisible(False)
    if item.text().replace('&','') == 'Database':
       item.setVisible(False)
    if item.text().replace('&','') == 'Web':
       item.setVisible(False)
    if item.text().replace('&','') == 'Processing':
       item.setVisible(False)
    if item.text().replace('&','') == 'Help':
       item.setVisible(False)

All menus (except 'Project' menu) are successfully removed. Only three menus that not succeed to be removed. They are 'Database', 'Web', and 'Processing'. Why are these menus not removed, while I have written the code for them?


texts: Enable saving only when specific columns were changed in QGIS 2.18

texts: I would like to force the user of a QGIS project file to edit 4 columns in the attribute table and only when those columns were edited the saving of the changes is enabled. Is this possible in QGIS 2.18?


texts: Are there any convenient sources of color palettes to use in ArcGIS software while producing maps?


texts: I am using QGIS 1.7.3 (can't get 1.8 to work on my WIN 7 64 OS) and I cannot find the Sextante plugin. I have found other plugins with no problems, accessing all the various repositories, I have a list of many to choose from several repositories but there is no Sextante although the Sextante website says it should be downloadable from the QGIS main repository. What am I missing, doing wrong?


texts: Apportioning Neighborhood Information to Census Tracts using ArcGIS Desktop?

texts:
I have a shapefile of census tracts and another shapefile of "neighborhoods" with values that I need to have joined to the census tracts. I know I need to apportion the area of the neighborhood shapefile to the census tracts (the neighborhoods are larger than the census tracts and I ultimately need the data by census tract) and then allocate the number of properties that paid property taxes in the neighborhoods to the census tracts proportionately. Hopefully I don't have that backwards.

Any suggestions?


texts: Developing ArcPad 10 forms with Visual Studio?

texts:
I have been trying to develop solutions with ArcPad 10 Studio for a while now and finding it restrictive.  I would like to use the Visual Studio environment to develop ArcPad 10 solutions (form development and coding).

Does anyone know if this is possible?


texts: Creating a fully custom map using OSM infrastructure

texts:
I would like to make a fully custom map (fictional, have to start from scratch basically) and use the OSM tools to make following things possible:

crowd editing - the source data which is used to generate the tiles (and in the end a large continuous map) will be publicly available, anyone can make a change to it (I don't require a nice web interface like OSM has, simple github repository will suffice)
custom rendering - I know that this is possible using custom tiles and such so this shouldn't be an issue.

Before I dive into this I would like to know if this is possible. I think tools like JOSM should be good enough for drawing the maps, then I can define my own tiles and use one of the rendering tools to display the map.

Are there any pitfalls I'm missing? Or maybe a better way to achieve this?


texts: Determining area amount of soil that needs to be excavated to level an area?

texts: I have contour lines, a raster, and a shapefile of the area that needs to be leveled. I have clipped the contour lines and the raster image to the shapefile of the area that needs to be excavated. I was thinking of creating a DEM from the contour lines, but I am not totally sure what to do next.


texts:
I am going to implement shapefile uploading from a webform, and ingestion into PostGIS.

One problem I am having is that my project is using Gradle for dependency management.

I am unfamiliar with Gradle and have no idea how to find the right string to get the jars into my project.


texts: QGIS 3 Proportional symbol by expression feature and legend display

texts:
When I use the tool's assistant in order to create proportional circles, I have no problem to display its legend. We can see the legend in canvas map and therefore, it's perfect for the print composer.

However if I use the expression feature, I have no possibility to display correctly the legend. The scale is not good. You can see the figure that values in legend are completely overestimated. It seems to be a unit problem or something else.

With the tool's assistant

With the expression feature


texts: How to fix undershoots in data converted from Microstation using ArcMap 10.0?

texts:
I have converted data from Microstation to ArcGIS. The converted data has undershoots at intersections, so I tried to snap lines using Spatial Adjusment. I've added some control points but ArcMap didn't fix them.

How can I fix the undershoots automatically? I have many lines in my data, and fixing them manually would be too time consuming.


texts: Validating topology using python script when ArcMap is open?

texts:
I am testing on validating topology using python.The situation is when I run the script when arcmap is open, an error occurs which indicates the topology has already been opened (FDO error: -2147214980 Cannot acquire a schema lock because of an existing lock), but when I close arcmap,the error disappear and it works properly.

Does it means I can't validate a topology when arcmap is open?


texts: Why is this ExtJS combo box not filtering?

texts:
I am trying to place a combo-box in the top of my web application which will search an internal WFS from Geoserver. I am using the combo-box so I can have the "type ahead" functionality which I think will help then narrow down the site they are looking for. Then once I solve this issue I will have a listener which will get the info and zoom to the location on selection.

I am working of a test at then moment to see if I can implement it in an existing site, but I am running into some trouble. The combo-box doesn't filter the data store, it just returns a full list of the datastore regardless of what i type (see picture at bottom of post).

The code I have used is below, unfortunately I don't have an external facing site to show.

I originally was hoping to set my WFS output format to JSON and just load that into a JSON data store but I kept getting errors that I was unable to decipher.

Thanks

Ando

var queryLayer = "dam";
var protocol= OpenLayers.Protocol.WFS({
    url: "/geoserver/ows",
    featureType: queryLayer,
    featureNS: "www.awebsite.com.au", //Just for SE proper featureNS website used in application
    version: "1.0.0",
    srsName: "EPSG:4283"
});

var assetStore = new GeoExt.data.FeatureStore({
    fields: [
        {name: "name", type: "string"}
        //{name: "type", type: "string"}
    ],
    proxy: new GeoExt.data.ProtocolProxy({
            protocol: protocol
        })
     });

new Ext.Viewport({
    layout: "border",
    items: [{
        region: "center",
        id: "mappanel",
        //title: "Map",
        xtype: "gx_mappanel",
        map: map,
        tbar:['->',{
                xtype: 'combo',
                fieldLabel: 'Search for an asset:',
                store: assetStore,
                displayField: 'name',
                forceSelection: true,
                triggerAction: 'all',
                minChars: 1,
                displayField: 'name'
        }],
        layers: [layer],
        extent: extent,
        split: true
    }]
});

Edit This is the query that is getting sent to the server, its not passing the name from the combo box to the query so I assume that why everything is being returned. Is there a configOption I am missing here??

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.0.0/WFS-transaction.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"><wfs:Query typeName="feature:dam" xmlns:feature="www.mywebsite.com.au"/></wfs:GetFeature>


texts: How I can use TSP function in pgRouting 2.0?

texts: I've spent years working with pgRouting and so far he had solved the problems I had needed. With the change to version 2.0, things improved, everything was much more organized and new algorithms that, which gives solution to dissimilar problems were incorporated. My question is related to the new TSP algorithm implementation, and not how to use it. I wish someone would put me an example of how on a network of 200,000 or more lines calculate a route with TSP of 5 0 6 waypoints.


texts: Finding point closest to point on near feature?

texts:
I have a dataset osmMajor containing road data (from OSM; 670,780 features) and another dataset points containing some points of interest (500,000 points). For each feature in my point layer, I want to find the coordinate of the closest point located on the road network.

I have two options to do this: either ArcMap's near tool (Analysis tools/Proximity/Near), or Spatialite. I tried using ArcMap, but it seems a very long task to perform (I had to stop the process after 3 hours). For this reason, I was thinking about using Spatialite to do this, but it seems that the process is less straightforward than what I thought.

How should I write my Spatialite query to get what I want?

For Spatialite, I think that I should follow these steps:

For each point feature, calculate the distance between the point and
the closest road
Around each point, create a buffer using the distance to the closest feature
Get the intersection between the buffer and the road data (it should be a point)
Get the coordinates (x,y) of the resulting point

Trying to implement these 4 steps on a sample dataset with only 10 points, and I encountered several issues.

For point (1.), I can find the overall smallest distance between the points and the roads, but I cannot have it for each point:

SELECT test.ID, min(Distance(test.Geometry, osmMajor.Geometry)) AS distance
FROM test, osmMajor
WHERE osmMajor.ROWID IN (
SELECT ROWID
FROM SpatialIndex
WHERE f_table_name = 'osmMajor'
AND search_frame = test.Geometry)
GROUP BY test.ID

I get a distance for only three points instead of the 10 that are in the test table. I suspect it has something to do with the way the spatial index works, because the query only returns the results for points 2, 3 and 4, which are the closest to the road network.

Can anyone give me advice on the best way to get the coordinate of the point on a feature which is closest to another point?


texts:
I got this error message when I tried to perfoming a spatial join.

An error has occured while executing Python code :
Traceback (most recent call last):
  File "C:/Users/Laurent/.qgis2/python/plugins\spatialJoin\spatialjoin.py", line 186, in applyJoin
    self.tra.ce( exp)
  File "C:/Users/Laurent/.qgis2/python/plugins\spatialJoin\spatialjoin.py", line 45, in ce
    print string
UnicodeEncodeError: 'ascii' codec can't encode character u'\xe7' in position 17: ordinal not in range(128)

Python Version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS Version:
2.8.1-Wien Wien, exported

Any idea on this?


texts: Where is Cluster point style in new Carto?

texts:
In the old interface for CartoDB (Carto Editor?) overlapping points could be displayed as clusters allowing users to click on a cluster to distribute the points visually and select one to read the pop-up. The Cluster style seems to be missing form the new Carto Builder interface.

Am I missing the style option for this or has this been removed from the new version of Carto?


texts: Converting models from QGIS to Python scripts?

texts:
I want to automate some tools for QGIS 2.0.2 using Python.

Previously I have used ArcMap for automating tools and now i want to use QGIS.

In ArcMap I used the ModelBuilder, and more specifically the command Export to Python Script to generate Python code from my model.

Now I want to find something similar that will do the same job for my QGIS model
i.e. I want to export my model in QGIS to a Python script.

Any ideas?


texts: For GPS point data using android, default camera app or other software?

texts:
I'm curious to know about your experience using the default android camera app for geotagging photos vs other apps.

I am going to trial a few different of these apps over the next month or so on different platfirms and I will share a report afterwards however I'm curious to know if anyone has any experience using the default camera app for point mapping and why you might recommend other apps. Do they have the advantage of increased spatial accuracy?


texts: Calculate area of all polygons in a shapefile and save it in attribute column "Area" in a Python script using QGIS

texts:
I have a shapefile with about 2000 polygons. I am using QGIS and I need to write a standalone Python script which opens the shapefile and calculates the area of all polygons separately and outputs the calculated area in a column "Area" in square kilometers. I tried to do it using the example in Calculate area in square meters from degrees?

eligibleLayer = QgsVectorLayer("/home/usr/Desktop/eligible_areas.shp", "eligible_areas", "ogr")

if not eligibleLayer.isValid():
    print "Eligible areas layer failed to load!"
lProvider = eligibleLayer.dataProvider()
lProvider.addAttributes( [ QgsField("Area",QVariant.Double) ] )
area = 0
for gFeat in eligibleLayer.getFeatures():
    calculator = QgsDistanceArea()
    calculator.setEllipsoid('WGS84')
    calculator.setEllipsoidalMode(True)
    geom = gFeat.geometry()
    landArea = gFeat.attributes()[gProvider.fieldNameIndex('Area')].toDouble()[0]
    if geom.isMultipart() is False: # if only simple polygon, calculate only for this
        polyg = geom.asPolygon() #transform to list of points
        if len(polyg) > 0:
            area = calculator.measurePolygon(polyg[0])

    else: #is Multipart
        multi = geom.asMultiPolygon()
        for polyg in multi:
            area = area + calculator.measurePolygon(polyg[0])

The script calculates the variable area but I am unsure if it is doing the right thing since I checked with the geometry tools and I get something completely different. Furthermore, I am unsure about how to update the above code such that I save the information in a newly created attribute column called "Area".


texts: Convert shapefile with lists to data frame in R

texts:
I would like to ask some help for converting a shapefile to data frame. I downloaded data related to roads from Natural Earth. I imported the shapefile using shapefiles package like the following.

foo <- read.shp("ne_10m_roads_north_america.shp")

Then,I checked str(foo) and saw the following. This is just the beginning part.

> str(foo)
List of 2
 $ shp   :List of 49183
  ..$ :List of 8
  .. ..$ record        : int 1
  .. ..$ content.length: int 656
  .. ..$ shape.type    : int 3
  .. ..$ box           : Named num [1:4] -119.6 37.8 -119.6 37.9
  .. .. ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
  .. ..$ num.parts     : int 1
  .. ..$ num.points    : int 79
  .. ..$ parts         : int 0
  .. ..$ points        :'data.frame':   79 obs. of  2 variables:
  .. .. ..$ X: num [1:79] -120 -120 -120 -120 -120 ...
  .. .. ..$ Y: num [1:79] 37.9 37.9 37.9 37.9 37.9 ...
  ..$ :List of 8
  .. ..$ record        : int 2
  .. ..$ content.length: int 40
  .. ..$ shape.type    : int 3
  .. ..$ box           : Named num [1:4] -119.8 39.1 -119.7 39.1
  .. .. ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
  .. ..$ num.parts     : int 1
  .. ..$ num.points    : int 2
  .. ..$ parts         : int 0
  .. ..$ points        :'data.frame':   2 obs. of  2 variables:
  .. .. ..$ X: num [1:2] -120 -120
  .. .. ..$ Y: num [1:2] 39.1 39.1

I usually do something like this.

df <- data.frame(matrix(unlist(foo), nrow=100, byrow=T))

But, it seems that each list had different numbers of rows. So I am not sure if this approach is right. I have also tried the following. do.call(rbind) will be probably slow even if it works. rbindlist() indicates the numbers of columns different among items. tidyr got a nice new function called unnest() which would work, but R crashed.

> test<- do.call(rbind.data.frame, foo)

Error in data.frame(record = 1L, content.length = 656L, shape.type = 3L,  :
  arguments imply differing number of rows: 1, 4, 79

> test <- rbindlist(foo)
Error in rbindlist(foo) :
  Item 2 has 12 columns, inconsistent with item 1 which has 49183 columns

It seems to be a relatively easy thing to do, but I am not familiar with shapefiles. How can I approach this situation?


texts: Looking for a tool for Orthorectification of Archived Air Photos

texts:
I was wondering if anyone knows of an open source, perferabily something for QGIS and not GRASS, that can be used for orthorectifying scanned archived air photos.  I have looked at OSSIM and a few other packages, but they do not allow (from what I can tell) for a single image to be processed based on user input from a geoimage or vectors.  Packages like ERDAS and PCI allow for this, but I am looking for an open source solution.  The Georefrencer plugin to QGIS is good, but I would like to be able to use a DEM (digital elevation model) as a source for the Z information for terrain relief issues.

Any help would be great, or clarification on OSSIM or other products that I may be misunderstanding.


texts: How to search within a string to extract data from OSM other_tags?

texts:
I am working with OSM data in QGIS and I often need to query items in the other tags field. When I download the data using the OSM  download plugin the other tags come in one field so are not easy to use.

Currently I am looking at permissive paths etc

"foot"=>"permissive","horse"=>"permissive","tracktype"=>"grade2"

I can use the expression (badly written) below to pull out what I want and perform a selection but it is a headache. any better ideas? I have tried quick OSM but that doesn't seem to have access to the details in the other tags

substr(  substr(  "other_tags" , strpos( "other_tags",'foot"=>') ,70),10,strpos( substr(  substr(  "other_tags" , strpos( "other_tags",'foot"=>') ,70),10,20),'","')) = 'permissive'


texts:
I'm running an update cursor on an in_memory feature class to remove blank (but not null) values as part of a larger process.

I wanted by update cursor to loop through my list of fields and if the value was '' or ' ' just replace it with a Null, however while it is finding the correct fields it's not updating the value.

If I explicitly refer to each item in the row e.g. row[1] it works as expected, but not if I try to loop through row.

This works:

with arcpy.da.UpdateCursor(AddressLayer, fieldList) as cursor:
    for row in cursor:
        arcpy.AddMessage(row)
        u = 0
        if i > 1000:
            break

        ##  THIS BIT DOESN'T WORK! ##
        #for r in row:
        #   if r == '' or r == ' ':
        #       arcpy.AddMessage("OID = {}, value = '{}'".format(row[0], r))
        #       r = None
        #       u = 1
        ## THE ABOVE DOESN'T WORK! ##

        ## ADDED BECAUSE THE ABOVE DOESN'T WORK ##
        r = 1
        while r < 9:
            if row[r] == '' or row[r] == ' ':
                arcpy.AddMessage("OID = {}, value = '{}'".format(row[0], row[r]))
                row[r] = None
                u = 1
            r+=1
        ## END OF ADDED BIT ##

        if u == 1:
            arcpy.AddMessage("Updating")
            cursor.updateRow(row)
        i+=1

But this (the commented section above) doesn't:

for r in row:
    if r == '' or r == ' ':
        arcpy.AddMessage("OID = {}, value = '{}'".format(row[0], r))
        r = None
        u = 1

What am I missing about the second code snippet that means it wouldn't work?


texts:
I am using OpenLayers to create a map where I can Zoom and export the view (Following this exemple https://openlayers.org/en/latest/examples/export-map.html)
This works fine, but what I need is a georeferenced image (ecw, geotiff) that users can use later on QGIS...

Does anyone knows how can I do this ?


texts: Re-project raster from dismo::gmap without losing colors

texts:
I downloaded a 'Google Map' terrain map using dismo::gmap. The raster comes with a Mercator projection. I tried to re-project it to UTM but in doing so, all the colors are lost and the raster is not displayed the same way (add scale + coordinates). Could you explain me how to retain the original colors ?

library(dismo)
library(raster)

area  <- extent(678815, 709272, 5004393, 5028388)
proj4 <- "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

r <- raster()
extent(r) <- area
proj4string(r) <- proj4

gm  <- gmap(x = r, type = "terrain", scale = 1)
gm2 <- projectRaster(gm, crs = proj4)

plot(gm)
plot(gm2)


texts: Linking pdf which is stored on local server to web map made with qgis2web

texts:
I have a QGIS project that I have linked pdf's to shapefiles using a file directory field with a "attachment"widget. I have it set up so when I click on the shape in my project, the link pops up which will bring me to the relevant pdf. I tried to export this project using qgis2web, but where the link normally is, there is just a question mark.

Any way to have that link in a web map? My intent is to not put the web map online, but have it on the internal server for co-workers to use.


texts: Is Spatial clustering with PostGIS possible to generate N number of groups that have 10 elements?

texts:
I need to create a query using some sort of spatial clustering algorithm within PostGIS database. Currently, I am using pgAdmin version 4. More specifically, I need to create N number of groups that have 10 elements. I know the output will end up with groups that have less but the groups should not have more than the 10 elements.

Can the variable for 10 be parametrized?

I have read the following post since it seems to provide useful information to achieve what I need but my level of expertise with PostGIS is not much. To only thing that I do not see that is addressed here is that I need to be able to select the number of groups and features to assign per group. It could be for point and/or polygon geometry.
Scenario:
I have a subset of 388 polygons. I need to make N number of groups with 10 polygons. The 10 polygons cannot be random but they must be adjacent or the closest possible. The output will be used to send inspectors to the field. Since they will be walking, it is preferred that the polygons shapes are adjacent to each other, or the closest to each other as possible.
I was able to achieve the following output using a QGIS plugin: ClusterPoints. It works great but it does not allow me to set the predefined number of features per group therefore, I ended up with groups with more than 10 features and others with three or four features assign to them.

Input:

Desired Output (this output has the flaws described above):

My intention is to have someone point me in the right direction.


texts: How to build cylinder or cube multipatch objects in ArcGIS?

texts:
I'm trying to visualise a series of boreholes, which consist of points on the ground surface, and soil sample values taken at 0.01m depth intervals (in a related table).

eg in this screenshot, I've selected one borehole and displayed its related samples in the table.

Ideally I'd like to create a series of cylinders in 1m depth segments, to represent the average soil values across that segment:

I can calculate the average soil values for any given 1m slice, and I know the start and end depth of each slice, and the desired width (eg 30cm).

My question: How can I build the cylinders?

I presume these will need to be multipatches, but I can't see anything in the Esri documentation on how to create the multipatches. Ideally I'd like to script this in arcpy as I have thousands of boreholes to process.


texts: Does anyone know a function to determine if a Point lies on a polyline?

texts:
ArcEngine 10 .Net 4  c#

I'm currently allowing user to select a point on a polyline and creating a point on the line.  I am using snapping tool.

I just want to make sure that when the user clicks on the screen, the point feedback that I get lies on top of the polyline.

I was wondering if there is a method that I can check if the IPOINT is ON any of the polylines in the polyline layer.


texts: How to get the intersections lat/lng on overpass api?

texts:
I've already knew that I can find all intersections on map by
Code.

My question is that I've knew one intersection and how I can get the information like lat and lon on others connecting intersections .

I've already get lat and lon of intersection A .
How to get the information on others connecting intersections in the following picture?


texts: PostGIS ST_MakeEnvelope size direction change on geography

texts:
This doesn't really effect me but I was just wondering if you find this a little bit strange or not. I think this issue is probably cause of this question also but I think the accepted answer does not apply to my question.

I am using ST_MakeEnvelope and it tax xmin,ymin,xmax,ymax values as arguments.

I realized that for example, if longitude distance is over 180 degrees, it seems to wrap around other way around (That is what I think is happening, I may be wrong?).

I thought it would always go from xmin to xmax. As in xmin=-81 to xmax=100 should go in direction -81,-80,-79...0,1...100 , yet it seems to make a boundary box going other way -81,-82,-83...(+/–)180,179...100

Am I imagining things and is this suppose to work like this? Any reason why this way is better than always assuming we should travel from xmin to xmin+ direction?

SELECT tbl1.column1, tbl2.column1, tbl1.column2 && tbl2.column2 AS overlaps FROM
( VALUES
(1, 'POINT(0 0)'::geography),
(2, 'POINT(-85 0)'::geography) ) AS tbl1,
( VALUES
(3, ST_MakeEnvelope(-80, -10, 100, 20, 4326)::geography),
(4, ST_MakeEnvelope(-81, -10, 100, 20, 4326)::geography)) AS tbl2;

Results in:

 column1 | column1 | overlaps
---------+---------+----------
       1 |       3 | t
       1 |       4 | f
       2 |       3 | f
       2 |       4 | t


texts: What is the difference in purpose between STAC and OGC services?

texts: I am confused about the purpose of the Spatiotemporal Asset Catalog. How does it differ from OGC services which specify accessing data through various services? Is it another service?


texts: Calculating area of user defined polygon in openstreetmap using python

texts:
My objective is to make a tool to calculate the area of any user defined polygon on the map using a python code. I tried doing this with OpenStreetMap.
My approach is:
Open a map and select points on the map which will be the vertices of the polygon. Export the lat long of these vertices into a python code which will calculate the area of the polygon. Integrate this python code into the map so that the area which is the output of the python code is displayed on the map.
Is this approach right?
The edit option allows us to draw polygons.

But when I try exporting this, it leaves the page. How do I get an xml file (or preferably the coordinates) of only the polygon that is marked?

I have marked a polygon in umap. How do I export the coordinates?


texts: Why Does My Raster Catalog Crashes When I Try to Zoom?

texts:
I have a file geodatabase which includes a Raster Catalog.  Up until today I've had no problems with my Raster Catalog.  Now, when I pull in the Raster Catalog, I can see the footprint of my rasters.  However, when I try to zoom to a specific area expecting the raster to show up, ArcMap crashes.

What's even weirder:  When I zoom into a specific area, the raster shows up.  But everywhere else, it crashes ArcMap.  What is going on??

I'm using ArcMap Desktop 10.2.2. The geodatabase is in version 9.3.  But even after upgrading it to 10, I still had the same issues. I've also tried the Recover Geodatabase Tool, but then it converts my Raster to a Polygon Feature Class based off of the Raster Footprints.


texts: Lines will not stop curving, can't draw a straight line

texts:
I am new to QGIS and I realize this is probably a really silly question but I am at my wits end I've been trying to solve this for hours. I need to make a simple line layer to represent the sewer system through town but my "straight lines" will not stop curving out!

What can I do about this?

I think it has something to do with the "add circular string" tool I just can't seem to turn it off?


texts: dynamic text added by default to every map document?

texts: I'd like to have the document path inserted by default on every map (that I can then delete if it's not needed) because I keep losing track of the electronic versions of maps for older ones (ca. 2012) I have printed, and my brain seems far too taxed to remember to add it each time. I have the standard ArcMap license with 10.3.


texts: QGIS atlas: data defined setting for individual scales is not working

texts:
I followed the following description of a data defined setting so that each page gets the predefined scale.

Atlas Scale Focus off

I have two "scale" coloumns with types integer and integer64 and the coverage layer is correct. However when I want to append the new scale field, I don't get any options to append.

What could be possible Errors? The setting makes perfect sense to me, it should work...


texts:
I used this to create line layer in the QGIS console:

iface.addVectorLayer("LineString", "pointlayer", "memory")

If I need a point layer what do i need to use?


texts:
Working in QGIS I have a point layer with ad attribute "NAME".
Is there any way to automatically create a folder called "NAME" where putting the attachments?
Maybe I could use something like this "How to check if a directory exists and create it if necessary?"
But how could I checking it directly from the active layer in QGIS?


texts: Clipping or masking a raster by individual polygons in a feature layer

texts:
I have a shapefile (polygon) layer consisted of four polygons and I have a raster. I want to calculate statistics for every raster part which is beneath any of these polygons.

I tried the following code:

import arcpy
import os
from arcpy.sa import *

os.chdir("E:/Reza/2015/spring and summer/1 oct/test")
inRaster="E:/Reza/2015/spring and summer/1 oct/test"+os.sep+"Band_1.tif"
objects="E:/Reza/2015/spring and summer/1 oct/test"+os.sep+"objects.shp"
with arcpy.da.SearchCursor(objects,["*"]) as cursor:
    for row in sorted(cursor):
        arcpy.CheckOutExtension("Spatial")
        outExtractByMask = ExtractByMask(inRaster, row[1])

but when I run, it says:

 Traceback (most recent call last):
      File "C:/Users/Rezars/Desktop/1", line 10, in <module>
        outExtractByMask = ExtractByMask(inRaster, row[1])
      File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy   \sa\Functions.py", line 7030, in ExtractByMask
        in_mask_data)
      File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy    \sa\Utils.py", line 47, in swapper
        result = wrapper(*args, **kwargs)
      File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 7026, in Wrapper
    out_raster)
      File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
      return lambda *args: val(*gp_fixargs(args, True))
     RuntimeError: Object: **[![Error in executing tool][1]][1]**

The table of my feature layer with its polygons is like this:


texts:
I'm using Copy Raster (Data Management) function. I need to ignore background values which are = 0 but keep all pixels =0 that are inside the raster.

background_value should do it, however it does not do anything. I found this explanation:

For file-based rasters, Ignore Background Value must be set to the same value as NoData for the background value to be ignored. Enterprise and geodatabase rasters will work without this extra step.

But when I set Nodata_value = 0 then I get rid not only of background values but all 0 values inside a raster. I'm not sure how this should help.
Code sample:

arcpy.CopyRaster_management(ras, final_data_downloads_dir + os.sep
                                    + output_file_name, background_value=0,
                                    pixel_type='8_BIT_UNSIGNED')


texts: Distance from points to points in the same polygon

texts:
I have 2 point layers and 1 polygon layer. I want to find distance nearest between this 2 point layer, but they must be in the same polygons.

I used distance nearest hub, NNJoin, MMQgis but they just have take part in of 2 layer, there are no conditions attached.

Apart from these tools, are there any tools or python or any way to solve this problem?


texts:
Im working with QGIS 1.8.0 and I have been trying to upload a delimited text layer but keep getting the following error :

The layer file:///C:/Users/fernando/Desktop/ATLAS/TAMAULIPAS/Estaciones%20Meteorol%C3%B3gicas%20Tamaulipas%20SMN.csv?delimiter=%5B,%5D&delimiterType=regexp&xField=X&yField=Y is not a valid layer and can not be added to the map

I attach you a picture of the settings Im using when uploading the CSV file. Do you have any ideas why could it be wrong?

Thanks in advance.


texts: How to get changeset tag from .osm file into PostGIS enabled Postgres db

texts:
End goal: visualize geographic data from .osm file with QGIS by adding PostGIS layers which include changeset information.

I have used osm2pgsql successfully to populate PostGIS and visualize the layers in QGIS but cannot bring in the changeset field.

I have used the --extra attributes function to bring in fields like osm_uid and osm_user but --extra attributes apparently does not allow you to bring in the changeset field.

How can I bring in the changeset field with osm2pgsql? Is there some alternative route I should take instead of osm2pgsql?


texts: How to get landmarks data contents that are added and deleted between 2 landmarks layers shapefile

texts:
Now I have 2 layers of landmark,2019 version layer(Blue) and 2020 version layer(Purple).
I want to know landmarks that are added and deleted in year 2020 and I see the landmarks that are added and deleted in year 2020 by order the layers.

How can I get the landmarks data contents that are added and deleted?
Thank you in advance.


texts:
I have a table with around 2 million records in it.  I create a spatial index, using the defaults other than the bounding box.  I've been noticing that some queries are extremely fast, and some are extremely slow.  The determining factor appears to the size of the polygon used in the query.

On larger search areas, using WITH(INDEX(SIX_FT5)) slows down the query considerably (from 0 seconds, to 15+ seconds).  On smaller search areas, the exact opposite it true.

Here's some of the queries I'm testing with:

Fast:

SELECT TOP(1000) * FROM [FT5] WHERE (shape.STIntersects(geometry::STGeomFromText('POLYGON ((-133462.805381701 -668610.241000959, 2934415.68824241 -668610.241000959, 2934415.68824241 2200521.65831815, -133462.805381701 2200521.65831815, -133462.805381701 -668610.241000959))', 2264)) = 1)

Slow:

SELECT TOP(1000) * FROM [FT5] WITH(INDEX(SIX_FT5)) WHERE (shape.STIntersects(geometry::STGeomFromText('POLYGON ((-133462.805381701 -668610.241000959, 2934415.68824241 -668610.241000959, 2934415.68824241 2200521.65831815, -133462.805381701 2200521.65831815, -133462.805381701 -668610.241000959))', 2264)) = 1)

Anyone know what's going on here?


texts: How can singleband pseudocolor style's values be computed using PyQgis

texts:
In Qgis desktop application there is a button named "Classify" to compute the values in singleband pseudocolor style. Highlighted in yellow in the following image.

How can I get those values using python?


texts: Using default datum transformation in QGIS?

texts:
QGIS 3.4: My project uses Pseudo Mercator as the CRS (datum: WGS 1984). I now added a shapefile which uses a local coordinate system from Austria (datum: MGI). I did not define a datum transformation in the general QGIS settings, neither did i define a datum transformation in the project settings->CRS (it is empty).

I would expect that no transformation is used and that there is an offset between my map and the shapefile (like in ArcGIS), but the shapefile aligns correctly with my background map (which is Pseudo Mercator, like the project CRS).

I assume that there is a default transformation, even if I did not define anything, but where can I see the default datum transformation or how does QGIS deal with that?


texts: ArcPad Fill Northing, Easting, Elevation

texts:
How do you get ArcPad to fill out the fields for Northing, Easting and Elevation of a point that was created/measured using the GPS on a hand-held unit?

The fields for are part of the geodatabase but I'd like to have them filled in when a point is measured in the field using the hand-held unit.

I have worked with shapefiles in the past but our new equipment is all geodatabase using .axf files for ArcPad. I don't know how to customize an .axf file.

I am a beginner when it comes to scripting and working with these files.

Any information or pointers would be extremely helpful!


texts: ArcGIS 10, ModelBuilder: calculate field - how to insert year values?

texts:
I have a series of polygon data. From each shapefile, I want to select the polygons with some condition (bb = 3).

To new shapefiles, I added a new field: "year".
I want to fulfill the "year" field with:

shp1 -> 2006
shp2 -> 2007
shp3 -> 2008
...
shp11-> 2014

I used ModelBuilder to iterate through shapefiles, make selection and add a new "year" field.

However, how can I calculate the year field by 2006,2007,...2014 values in iteration mode?

Or, how to at least insert in year Field unique values? I tried to insert there my %name% variable (as mlc_2006, thus not the number but a text),

but even the "field" is as TEXT, I have an error:

I have no experiences with Python, thus I prefer ModelBuilder.


texts:
I am new to arcpy and I'm trying to write a script to get the raster and check for the minimum cell size value using Get Raster Properties method.

After I run my code I get a message of the minimum value but the raster is not saved. Below are my line of codes:

import arcpy ,os ,sys
from arcpy import env
from arcpy.sa import *

input_raster= arcpy.GetParameterAsText(0)
output_folder = arcpy.GetParameterAsText(1)

arcpy.CreateFileGDB_management(output_folder,"output.gdb")
outputgdb = output_folder +"\\output.gdb"+ os.sep

output_con_raster = outputgdb + "cond_raster"

inRas1 = arcpy.Raster(input_raster)
rasterValue = arcpy.GetRasterProperties_management(inRas1,'MINIMUM')
elevation_min= rasterValue.getOutput(0)
elevation_min.Save(output_con_raster)

Is there anything wrong with my code as the save method is not working?


texts: Hyperlinking logo in ArcGIS Web AppBuilder?

texts:
I would like to insert a hyperlink into my organizations logo, in the upper left corner of my ESRI Web AppBuilder app (see photo). Clicking on it currently opens a link to a blank web page. It is hosted on our web server, rather than ArcGIS Online.

Can anyone point me to the correct file and code snippet to edit it?


texts: Adding timestamp to NetCDF xarray dataset

texts:
I am new to working with NetCDF datasets and Xarray in Python. I have an xarray dataset made from 100 TIFF raster files:

As you can see there are 100 data variables, that are essentially the individual tifs:

Now, I want to add time as an additional dimension with band, x, and y. Since I have 100 variables, I created a time series with 100 timestamps (just as a random test):
timestamps = pd.date_range("2020-01-01", periods=100)

But when I try to add it as a new dimension, it throws and error RuntimeError: NetCDF: HDF error.
What do I need to do to get this to work?
merged.expand_dims(time=timestamps)


texts: Automating calculation of distance between zip code (centroid) and address within one line of data?

texts:
I have a large data set of federal contracts which includes two pieces of spatial data per entry (or line): the zip code in which the project took place and the address of the contractor.  I am interested in automating the calculation of distance between these two points for each contract and storing it in a new field in the attribute table so that I can run some simple and complex statistical analyses.

To further complicate the matter, I need to calculate distance in road miles and travel time (which i have done in other projects using arcgis online analysis tools).  I have never figured out how to do it in ArcMap.


texts: Openlayer3 basemap appears white block in Chrome

texts:
i faced interesting an issue on Chrome. im using google's 4 number base layer (street, hybrit, etc.) When my tiles are loading over baselayer, baseLayer becomes white block. Why?

No any issues on Firefox and IE8 or IE9.

 var map = new OpenLayers.Map('map', {
        projection: 'EPSG:900913',
        layers: [
            new OpenLayers.Layer.Google(
                "Google Physical",
                { type: google.maps.MapTypeId.TERRAIN }
            ),
            new OpenLayers.Layer.Google(
                "Google Streets",
                { numZoomLevels: 20 }
            ),
            new OpenLayers.Layer.Google(
                "Google Hybrid",
                { type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20 }
            ),
            new OpenLayers.Layer.Google(
                "Google Satellite",
                { type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 20 }
            )
        ],
        center: new OpenLayers.LonLat(35.5, 39.0)
            .transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913")),
        zoom: 7
    });


texts:
I installed QGIS 2.14 on a Windows7 64bit via the installer. But QGIS can not boot up remaining "Restoring loaded plugins...".

Although I re-installed it after I un-installed all of QGIS software including older versions and removed QGIS-related folders and registory key once from the pc, the 2.14 still can not boot up.

Does anyone know a solution?


texts:
The following script only gives me a result for the last given grid (grid2"), does anyone know why?:

# Import arcpy module
import arcpy
import shapefile
# Arcpy settings
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
# =======================
lstgrids = [[grid1],[grid2]
            ]
# loop
for i,item in enumerate(lstgrids):
    A,B = item

# Process: Project Raster
print 'project'
arcpy.ProjectRaster_management(A,B, "PROJCS['Belge_Lambert_1972',GEOGCS['GCS_Belge_1972',DATUM['D_Belge_1972',SPHEROID['International_1924',6378388.0,297.0]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Conformal_Conic'],PARAMETER['False_Easting',150000.01256],PARAMETER['False_Northing',5400088.4378],PARAMETER['Central_Meridian',4.367486666666666],PARAMETER['Standard_Parallel_1',49.8333339],PARAMETER['Standard_Parallel_2',51.16666723333333],PARAMETER['Latitude_Of_Origin',90.0],UNIT['Meter',1.0]]", "NEAREST", "100 100", "'Amersfoort_To_ETRS_1989_MB + Belge_1972_To_ETRS_1989_2'", "", "PROJCS['RD_New',GEOGCS['GCS_Amersfoort',DATUM['D_Amersfoort',SPHEROID['Bessel_1841',6377397.155,299.1528128]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Double_Stereographic'],PARAMETER['False_Easting',155000.0],PARAMETER['False_Northing',463000.0],PARAMETER['Central_Meridian',5.38763888888889],PARAMETER['Scale_Factor',0.9999079],PARAMETER['Latitude_Of_Origin',52.15616055555555],UNIT['Meter',1.0]]"


texts:
I still have trouble reprojecting shapefiles :(

I would like to reproject it to UTM33N / WGS84 (EPSG 32633) I think the projection of the shapefile might be Pulkovo 1984/ Gauss-Krüger Zone 4 (EPSG 2398) or ETRS89 / UTM33N (EPSG 25833) http://projektidee.org/images/extern/ffh.zip

I failed to reproject it though...

Could anyone here help me with that?

I'm trying to do this with QGIS.

What I do is right click on layer an save as - then I choose the Projection system UTM33N/WGS84 and save the file to a new shapefile.

The other thing was to first assign the other projection (pulkovo) to the shape layer and then save it to a new shapefile.

With both ways of dealing with it I didn't succeed...


texts: Finding vector data for UK railway lines?

texts:
Can someone point me in the direction as to where I can obtain vector layers for the following UK railway lines:

DLR,
GREAT ANGLIA,
TfL RAIL,
OVERGROUND,
C2C,
SOUTH EASTERN


texts: Disable saving vector layer in QGIS, especially in DXF and KML file?

texts: Can you help me on how to disable the exportation of vector layers especially in AutoCAD DXF, DBF File and ESRI Shapefile?


texts: How can I synchronize Google Maps, Wikipedia, Mapmyindia, and Bing Maps data on one webpage?

texts:
I want to combine maps from Google maps, Wikimapia, Mapmyindia and Bing maps on one web page with synchronization like in Dual Maps.

For example, if I pan on Google maps it should show the same location on all maps.

Does anyone know how to do this?

Any help is appreciated.


texts:
What parameters should be set (x spacing, y spacing, units) to get a 1 ha hexagon? I work in EPSG 32634 projection. QGIS plugin MMQGIS.

Regards,

Onosz


texts: Not found the code for a map R animation

texts:
I tried to make an animated map in R with tmap package, like here
https://geocompr.robinlovelace.net/adv-map.html#animated-maps

I tried to find the code for the animated map presented in section [8.3 Animated maps] from the link but I have not find it.
Where is the code for this map?


texts:
I use PROJ.4 cs2cs to convert data from a Lambert Conformal Conic to a geodetic coordinate (NAD83) or vice versa. In my application, sometimes I need to convert a lot of data. Because I call cs2cs with subprocess.Popen() and parse result in python, input/output overhead becomes too big for a large data set.

Is there an faster alternative to cs2cs or a way to make calling cs2cs in python faster?


texts: Distinguish between agricultural lands and other land types in the Google Earth Engine

texts: I have classified an image into 5 different categories of soil land, build-up, jungle, water, and the lands that are green but are not jungle. Now, I want to write a code to separate and distinguish the agricultural lands from the last category, so that I can classify my image to 6 categories. Does anyone know how to do that?


texts: How to download sentinel 1 for a particular polygon from Google Earth Engine

texts:
I want to download Sentinel-1 SAR data for both VV and VH with IW only within a rectangular polygon? I could download with the ROI, but even though ROI is a rectangle, the output to export is a rectangle.

The red rectangle is my ROI, but it's taking the extreme coordinates and making a square that encloses the ROI and exporting that output as shown.
I only want the data inside the ROI.
I even tried to use ee.Image.clip. But I couldn't able to comeup with a solution which can download the data. It was going to a large file
var start_date = ee.Date('2019-05-01');
finish_date = ee.Date('2019-06-15');
var orbit = 'ASCENDING';

var collectionS1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    // .filter(ee.Filter.eq('orbitProperties_pass', orbit))
    .filterDate(start_date, finish_date)
    //.filterBounds(polygons);

 // Get the VV collection.
var collectionVV = collectionS1.select('VV');
// Get the VH collection.
var collectionVH = collectionS1.select('VH');

var VV = ee.Image(collectionVV.first().clip(roi_A);  //roi_A is my rectangular roi
var VH = ee.Image(collectionVH.first());

// I would like to take VV and VH and stack them and download
var VVVH = VV.addBands(VH)

Export.image.toDrive({
  image: VV,
  description: 'A_ERT',
  scale: 10,
  //region: ROI,
  fileFormat: 'GeoTIFF',
  //formatOptions: {
  //  cloudOptimized: true
  //}
});

I'm getting the following error message.

Error: Export too large: specified 6875185131 pixels (max: 100000000). Specify higher maxPixels value if you intend to export a large area.


texts: Configuring output table popup for geoprocessing widget

texts:
I have built a simple geoprocessing widget that produces a table (record set) as its final product. I am trying to simply change how and which fields are displayed. I am using Flex Viewer 3.5 and Arc Server 10.1.

Here is my output:

Basically I want to only show the Value and Acres field. I tried doing this through the flex viewer in while configuring the GP Widget, but it didn't have any popup configuration options.

Any ideas of where I can set these configurations? GP.xml perhaps?

Thanks in advanced.

Adam


texts: How to add an image (jpg) to an existing layout using Python?

texts:
I have 1500 MXDs with the same template layout, and need to insert other images for the all MXDs, can I do this using Python?

I'm using ArcGIS 10.1


texts: Using Google API within Python and applying to ArcMap?

texts:
I am currently doing a practical project at university (MSc in GIS).

I have to develop a python script and the functionality of this script has been left very open other than we have to use Python only and with ArcMap.

I’m have use ArcMap and Python which I am very new to.
I want to create a popup tool\form that uses python – the tool\form would allow allow me to enter an Address in a text field.
When I press a Go button the python script uses GeoPy (Python Geocoding Toolbox)  to query google for the coordaintes of this address.
These coordinates and the address details are shown in the tool\form and the point is added to my active layer.

Is a Python Toolbox the correct functionality within ArcMap to use?
Is it possible to import the GeoPy scripts into the Python tool script using the ‘import GeoPy’ command (should I have the GeoPy scripts in the same folder as the Toolbox file)
What if multiple addresses are returned from Google – does the Python Toolbox allow me to view the array of address and select one – i.e. does the python toolbox allow user interaction or should I be using some other functionality  of ArcMap to do this with python?


texts:
I'm calling the following URL through a JavaScript application to obtain the unique neighborhood's names of my neighborhood layer.

http://.../geoserver/wps?service=WPS&version=1.0.0&request=Execute&identifier=gs:Unique&DataInputs=attribute%3dnome;features%3d%40xlink%3ahref%3dhttp://.../geoserver/smurb/ows%3Fservice=WFS%26version=1.0.0%26request=GetFeature%26typeName=smurb:bairros_vigentes_3857%26propertyName%3dnome&RawDataOutput=result%3dformat%40mimetype%3dapplication%2Fjson

The problem is that this call takes like 10 seconds to execute and the layer has like 90 geometries only. If I make the "same" call through GeoServer's WPS Request Builder, it executes almost instantaneously, but I suspect it doesn't do HTTP requests.

Anyway, there are any tips to improve it's performance? Other suggestion to obtain the same results are also accepted :)

Additional info:

Geoserver 2.7.2
Linux machine
"..." in the URL means my GeoServer host


texts:
I'm trying to export a series of images using the TimeManager plugin.

I'd like to have a legend on each, to demonstrate why my datapoints are sized the way that they are (i.e. this circle size represents one square foot, this circle size represents ten square feet).

Is it possible to add and style a legend in base QGIS? I do not want to use Print Composer as that won't work well with TimeManager.

I am fairly new to QGIS.


texts: How to use external CartoCSS/mss within JavaScript?

texts:
I am using CartoCSS to style features in a run-time CartoDB layer (as opposed to a pre-made visualization) but it gets messy doing this within JavaScript:

sublayers: [{
  cartocss: '#initvisitedfws {bottom/marker-width: 24; bottom/marker-ignore-placement:false; bottom/marker-fill: gold; bottom/marker-allow-overlap:true; top/marker-width: 14; top/marker-fill: #111; top/marker-file: url(clock.svg); top/marker-allow-overlap:true;',
  interactivity: 'cartodb_id, parent',
  sql: "SELECT * FROM subprojects"
}]

I could keep it somewhat neat by using multi-line strings but that becomes unsightly and tedious in its own way, so I would prefer to edit the CartoCSS in its own native .mss file (especially in tandem with this Sublime Text 2 .mss syntax/autocomplete package).

The ghetto-ish way I am approaching this is with an XMLHttpRequest:

function getMss() {

  var xmlhttp;
  xmlhttp=new XMLHttpRequest();
  xmlhttp.open('GET', "somepath/somefile.mss", false);
  xmlhttp.send();

  var mapStyle = xmlhttp.responseText;
  return mapStyle;

}

...and then calling it from the CartoDB sublayer parameters:

var cartodbParams = {
  user_name: 'user_name',
  type: 'cartodb',
  interaction: true,
  sublayers: [{
    cartocss: getMss(),
    interactivity: 'cartodb_id, parent',
    sql: "SELECT * FROM subprojects"
  }]
}

This works fine but leads to another HTTP request, and I'm also not optimistic about its maintainability once I have dozens or hundreds of rules and multiple map layers.

So, what is a better solution for using CartoCSS within JavaScript?


texts: What transformation is used for EPSG 28409 in GDAL?

texts: I used gdalwarp to convert a raster file in EPSG:4326 to one EPSG:28409 I believe that using this method doesn't apply a full transformation as there is no option to select between the numerous transformation options available for 28409 so the output is shifted from what I am expecting.  Does anyone know about this and if there isn't a way to do it using gdalwarp is there a way to do this using another gdal tool?


texts: Count number of different occurrences in a cell

texts: I have points that correspond to the observation of different species. These point are grouped in taxonomic orders (for those not familiar with the term, just imagine that a few hundred species are aggregated into a dozen or so groups). I want to make a raster in which the value of the cells corresponds to the number of different orders in each cell. I've tried spatial join and a few different things, but I haven't managed thus far.


texts: Gdal image georeferencing and gdal2tiles image flip question

texts:
I'm georeferencing an image using gdal and the following command:

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0.016264915466308594 51.505109712517786 -0.0024890899658203125 51.492018739579336 -a_srs EPSG:4326 Penguins.jpg penguins.tif

then I'm generating the tiles using gdal2tiles and the following command:

python gdal2tiles.py -t 'Penguins' -k penguins.tif penguins_gdal

but for some reason the resulting tiles show the image flipped on the y axis even though the tif image is not flipped.

Am I doing something wrong?

Thanks


texts: Join field repeatedly on multiple tables

texts:
I have a fifty or so shapefiles. All of them have similar attribute tables with identical field titles. I also have a standalone table with data that I want to add to every single attribute table based on a shared field (DATE).

This is easy individually using the "Join Field" tool, but on my computer it takes a long time (almost seven minutes per file), and I have to stay near the computer to start the next file.

Is it possible to set it up to run on every single shapefile in the batch? I want to let it run while I'm away from the computer.


texts:
Im working with QGIS 1.8.0 and wanted to ask you how can I assign to each of the small squares of the map (I have the coordinates of the center each of them) values from an interpolated data layer as the one shown below? I tried to use the point sampling tool, but I can not save the interpolated data layer (at least not with the right mouse click option) so it doesn't appear in the "Layers with fields/bands to get values from" text box.

Any ideas?

Thanks in advance!


texts: Custom DescribeFeatureType response in Geoserver?

texts:
I'm using Geoserver in my web application. I need custom DescribeFeatureType response. What this mean?
layer_a is a layer in Geoserver. It have Feature Type Details as follow:

DiscribeFeatureType request for this layer, return a response as follow:

<?xml version="1.0" encoding="UTF-8"?>
  <xsd:schema xmlns:gml="http://www.opengis.net/gml"
              xmlns:opengeo="http://opengeo.org"
              xmlns:xsd="http://www.w3.org/2001/XMLSchema"
              elementFormDefault="qualified"
              targetNamespace="http://opengeo.org">
  <xsd:import namespace="http://www.opengis.net/gml"
              schemaLocation="http://localhost:8080/geoserver/schemas/gml/3.1.1/base/gml.xsd"/>

  <xsd:complexType name="WorldCountriesType">
    <xsd:complexContent>
      <xsd:extension base="gml:AbstractFeatureType">
        <xsd:sequence>
          <xsd:element maxOccurs="1" minOccurs="0" name="the_geom" nillable="true" type="gml:MultiSurfacePropertyType"/>
          <xsd:element maxOccurs="1" minOccurs="0" name="ObjectID" nillable="true" type="xsd:int"/>
          <xsd:element maxOccurs="1" minOccurs="0" name="FIPS_CNTRY" nillable="true" type="xsd:string"/>
          <xsd:element maxOccurs="1" minOccurs="0" name="SQKM" nillable="true" type="xsd:double"/>
          <xsd:element maxOccurs="1" minOccurs="0" name="SQMI" nillable="true" type="xsd:double"/>
          <xsd:element maxOccurs="1" minOccurs="0" name="COLORMAP" nillable="true" type="xsd:int"/>
        </xsd:sequence>
      </xsd:extension>
    </xsd:complexContent>
  </xsd:complexType>
  <xsd:element name="WorldCountries" substitutionGroup="gml:_Feature" type="opengeo:WorldCountriesType"/>
</xsd:schema>

I want to get more description for some fields. For example, for SQKM row, I want to get some extra description such as: minValue, maxValue and etc.

Is it possible? If yes, HOW?

A suggestion: We can put our data in column comment of table and then define a db datastore in geoserver. So we can response comment of table with describeFeatureType?


texts:
I'm fairly new to GIS.

I'm using QGIS to overlay a vector grid on my map and getting data on different attributes within each box of the grid. I need to do this for different grid sizes hence there will be some repetition. So to make this efficient I wanted to write a PyQGIS script to run a loop to do this. However, I'm new PyQGIS too and am not familiar with the commands (or the library itself).

Is there a way to get the PyQGIS code when I implement something from the QGIS GUI?

In case someone is familiar with STATA, I'm looking for something analogous to getting the STATA code when you use the point-and-click method to run a command.


texts: I am using mapfish module to print my map to pdf. It works well for two layers but if I add another layer there is no response from geoserver. Neither there is any error in geoserver log no in apache.


texts: Postgis Geo Elevation Raster: Get average Elevation

texts:
I am running Postgresql 9.5.5 with Postgis 2.2.1 on Ubuntu. I have imported ASTER GDEM elevation information as a raster and query it for elevation using geo coordinates.
Currently I use

ST_NearestValue(elevation.rast, ST_SetSRID($POINT,4326), false)

With (e.g.)

$POINT = ST_SetSRID(ST_MakePoint(8.4255, 49.0113), 4326)

but as ASTER GDEM is "only" sampled at 30m x 30m I need to get a (weighted) average of the nearest 4 or 16 available samples. Which is the most elegant (and performant) way to do this?

I tried

SELECT (ST_Intersection(elevation.rast, ST_Expand($POINT),0.0001))).val
FROM elevation
WHERE ST_Intersects(elevation.rast, ST_Expand($POINT),0.0001));

but this only returns the height without the exact position, so the samples can not be weighted by distance. It also seems rather ugly and the expansion factor needs to be calculated depending on the latitude.

I am rather new to postgis and may have missed even a simple solution.


texts: ArcGIS Portal - Click to add a marker and get coordinates

texts:
How would this be best approached in ArcGIS Online and/or Portal?

A user needs to click a point on a web map to identify a work location.

Clicking the map adds a marker and also stores the XY coordinates in a variable to send to a external apps, databases and files.


texts: Assure the continuation of rivers (raster dataset)

texts:
I am creating a cost surface raster for identifying least cost paths for road constructions. One of the factors influencing the cost is water cover. I am using Landsat based Global Surface Water raster dataset for water cover. In some cases the rivers on the raster dataset do not have continuation and they break, messing up my least cost paths (assuming that on the break point there is no water, but there is).
How can I assure the continuation of rivers?
I am using ArcGIS Desktop 10.3.