texts: How to create a raster dataset from a known decimal extent?

texts:
I have an input raster pInputRaster. I want to create a new raster dataset called pOutputRasterDS that is the same width, height and cellsize as my input raster dataset (I'm not worried about populating the new raster with values at this moment in time). The input raster dataset has square cells.

I can get the input raster dataset information like so:

Type factoryTypefinalraster = Type.GetTypeFromProgID(
                            "esriDataSourcesRaster.RasterWorkspaceFactory");
IWorkspaceFactory pWSF = (IWorkspaceFactory)Activator.CreateInstance //DECLARE RASTER WORKSPACE FACTORY. THIS IS A SINGLETON OBJECT. ONLY ONE OF THESE PER PROGRAM
                            (factoryTypefinalraster);

ESRI.ArcGIS.Geodatabase.IRasterDataset2 pInputRasterDS = (IRasterDataset2)pInputRaster

ISpatialReference sr = new PointClass();
ESRI.ArcGIS.Geometry.IEnvelope pREnv = pInputRaster.Extent;
double dLowerLeftX = pREnv.LowerLeft.X;
double dLowerLeftY = pREnv.LowerLeft.Y;
origin.PutCoords(dLowerLeftX, dLowerLeftY);
int NumBand = 1;

ESRI.ArcGIS.Geodatabase.IRasterWorkspace2 pTempRWS2 = (IRasterWorkspace2)pWSF.OpenFromFile(strOutRasPath, 0);
string outfilename = "outputRas";

Now I have all that I need to call IRasterWorkspace2.CreateRasterDataset(). Here I call that method:

IRasterDataset2 pOutputRasterDS = (IRasterDataset2)pTempRWS2.CreateRasterDataset(outfilename, "GRID", origin, pREnv.Width, pREnv.Height, cellsize, cellsize, NumBand, rstPixelType.PT_UCHAR, sr, true);

The issue arises in the implementation of this method. The method takes type int or type long as inputs for the width and the height. In my case, I am most often working with a geographic coordinate system with decimal degrees, so my width and height are of type double. Obviously, converting these values to int is not possible. Projecting the data beforehand is also not a possibility. Any alternate ideas? I have considered just making a copy of my current raster dataset and clearing all values of it, but have not explored this option further.


texts: QGIS Regression Tools using only feature data

texts: I'm trying to find in QGIS an analogous tool to ArcGIS's "Ordinary Least Squares".  The tool takes in a shapefile, where the dependent variable and independent variables are all fields in the same shapefile.  The tool runs OLS and returns a shapefile with the same features as the original, with fields for fitted values and residuals.  The SAGA tools in QGIS are not very intuitive nor well-documented, but from what I can decipher none of them seems to work like the "Ordinary Least Squares" tool in ArcGIS.  To be clear, I want to run a regression where none of the data is in grid/raster format.


texts: Google Earth Engine, Creating a Feature Collection (Geometric points)

texts:
I have created an image from the crop landcover dataset containing only the crop type rice (crop_rice). I hoping to make a feature collection (Geometric Points) on all points containing rice. Any ideas on how to go on doing that?

//Select first image crop landcover data set
var dataset = ee.ImageCollection('USDA/NASS/CDL')
              .filter(ee.Filter.date('2017-01-01', '2018-12-31'))
              .first();

//select band cropland from image
var img_cropLC = dataset.select('cropland');
Map.setCenter(-100.55, 40.71, 4);

//Add visualization layer to view all crop types in the image
Map.addLayer(img_cropLC, {}, 'Crop Landcover');

print(img_cropLC)

//Select all pixel with crop type:Rice (crop class values = 3)
var crop_rice = ee.Image(img_cropLC).eq(3)
print(crop_rice)

//Add a visualization layer to view rice
Map.addLayer(crop_rice.mask(crop_rice),{
bands:'cropland',
palette:'0000FF',
},'Rice Layer')


texts: Converting MULTIPOLYGON field stored as text to geometry data type in PostGIS and plotting in QGIS


texts:
I am working on a tool, which transfers data, including geometries, from different legacy databases to postgis.
Currently this means I have a MS-access Database, which contains point geometries, but since access does not have a geometry type, the points consist of two columns, containing numbers, like so:

lat
lon

371744,668
5700069,115

During the flow of my tool, I take these coordinates with pyodbc, and use shapely to create a point. This point does not have crs since it is newly created. Now I'm trying to set it to EPSG:25832, since the postgis DB expects geometry data to be in 25832. The code below shows what is going on. In the docs I found an example, and even when I'm directly pasting the example to my code it does not work.
What am I missing here?
def getGeomAsWKB(self, table, ids):
    tmp = self.getCoordinates(table) # tmp = list of tuples with coordinate pairs [(X,Y),(X2,Y2)...]
    geoms = []
    for coords in tmp:
        pt = Point(coords) no crs at this point, as to be expected
        srid = pyproj.CRS('EPSG:25832')
        old_srid = pyproj.CRS('EPSG:25832')
        project = pyproj.Transformer.from_crs(old_srid, srid, always_xy=True).transform
        t_pt = transform(project, pt) nothing changes
        #pt._crs = 'EPSG:25832'# i tried to set it this way, becausei did not know what else to try
        tmp_binary = wkb.dumps(t_pt)
        geoms.append(tmp_binary)
    return geoms

Just to be clear, pt and t_pt are correctly build, its just that t_pt does not contain the crs information after the transformation.
Postgis return the following error code :
Geometry SRID (0) does not match column SRID (25832)


texts: Seeking Trimble GGF Geoid File Format Specification?

texts:
I want to generate Trimble GGF Geoid files from my software.

Can anyone provide me with a format specification for the binary ggf files used by Trimble?

I found a basic description on the OpenStreetMap wiki:

ORIGIN: top left
0 0 TNL GRID FILE\0 XXX\0*
1 0 TNL GRID FILE
0x30 8-> ymin (°) (Min Latitude in decimal degrees)
0x38 8-> ymax (°) (Max Latitude in decimal degrees)
0x40 8-> xmin (°) (Min Longitude in decimal degrees)
0x48 8-> xmax (°) (Max Longitude in decimal degrees)
0x50 8-> dy (°)    (Step size for latitudes, in decimal degrees)
0x58 8-> dx (°)  (Step size for longitudes, in decimal degrees)
0x60 4->  (vNY, nrows)
0x64 4->  (vNX, ncols)
       North Pole
       South Pole
       Missing Value (nodata)
       Scale
       Flags*
0x92+ data*

But some information is missing:

How is the data between 0x68 and 0x92 supposed to look?
How are the geoid heigth's supposes to be stored in the data* section?

Column-wise?
Row-wise?
South-West to North-East or the other way around?


texts: Boolean fields from postgis data type changed to QString in QQIS

texts: I have a layer in QGIS loaded from a postgis database. In the layer properties dialog in QGIS, under the fields tab, layers with a typename of 'bool' (which are boolean fields in the postgis database) appear as type 'QString'. The values then behave as strings in QGIS. Is there a way to bring the layer in so that the boolean fields are treated as booleans in QGIS?


texts: Contouring with cliff faces in ArcGIS for Desktop?

texts:
I am trying to create some contour lines that are around cliff faces.  What is the best way to create contour lines that are effected by cliff faces and is there a way to automatically draw the cliff faces?

I've tried using the Contour with barriers tool in ArcGIS but can't seem to get it to work.  I've downloaded the DEM Raster and can produce contours from that, just can't take it to the next step where the contours are effected by barriers.

I'm using both ArcGIS and QGIS.


texts:
We have spatial layers in ArcSDE with objectid as unique value, I've taken this data to another platform & migrated data to cloud.

Now, after months, back in ArcSDE and want to migrate new data (huge data), but want to take only DELTA data to cloud, i.e. any new records, deleted records and modified records, so I can insert new records into existing data in cloud, delete records which have been deleted in source, and update records which are modified in source.

So, my questions are:

Is OBJECTID in feature class in SDE are reliable to calculate whether they are new, modified & deleted? Can I compare earlier export of the data & find objectid's?
What are options in SDE or outside (examaple: PostGres) to check only delta records?
Do I need to separately enable in SDE to capture delta or create additional fields to capture D,I & U records?  Our geodatabase is versioned, and multiple people are creating versions and being posted to production regularly.


texts: How to load a raster from in-memory bitmap (not from file) into PostgreSQL

texts:
I have a bitmap in memory that I would like to load as a raster to PostgreSQL.

An unfavorable option is saving the bitmap to file and then use raster2pgsql, but as we
have several cameras that produce images at high rate, we can't afford saving them to files just for sake of insertion to the DB.

I'm looking for a solution for loading the bitmap directly from memory. Any suggestions?

I'm using PostgrSQL 9.2 using PostGIS2.0 on Windows 7.
The images are plain gray level, if it matters.


texts: Legend of Flow Direction calculated with processing toolbox QGIS (SAGA)

texts: I tried to figure out the legend of the output for flow direction when I use the QGIS processing toolbox function from SAGA (Channel network and drainage basin). I get a raster with values from 0 to 7. The Taudem d8 uses an encoding from 1 to 8 that does not correspond with this result. Anybody knows how this is encoded?


texts: Comparing TIGERLine 2006 (and prior) With New Shapefile Format

texts:
I have created a geo database using MySQL and the TIGERLine data files available from the US Census.  This is actually an MSSQL database I created that I've transfered to MySQL.  Everything works beautifully for geocoding and reverse geocoding.  The only problem is that data is now almost 6 years old and I am looking for greater accuracy.

I've looked at the 2010 TIGERLine files and noticed that they changed the format back in 2007 to a new Shapefile format. I've had some success using Ogr2Ogr in getting those shapefiles populated into a database but now I'm stuck.  I was able to build my street segments table (center line data) easily with the 2006 data and the format made sense from a processing point of view but the new shapefile format doesn't seem as straightforward.

Is there any documentation (and I've scoured the US Census Bureau) that relates the two formats?  A side by side comparison of fields would be too much to ask, I'm sure, but that is exactly the kind of data that I'm looking for.

Thanks in advance.


texts: Thematics in Mapinfo - multiple intervals from same table?

texts: I´ve got a table where I want to thematize intervals of values like between 0.1 - 1.0, 1.0 - 2.0 and etc., but I´m almost sure that´s not possible, whithout splitting up the table to multiple tables, right? It has to be for Mapinfo.


texts: Looking to do a one-side buffer with QGIS v2.18.9

texts: I'm working on a project that requires a 30m buffer around a linear feature. The feature represents the coast and I need the buffer to show 30m inland from the line representing the coast. Is there a way to do this using QGIS?


texts:
I have a date field in attribute table looks like this "11/04/2020 14:25:05" and I am trying to run the below function from field calculator

def date_info(strDT):
    import datetime
    dt_obj = datetime.datetime.strptime(strDT, '%d/%m/%Y %I:%M:%S')
    weekday = dt_obj.isoweekday()
    return weekday

but it returns an error as below
any ideas?


texts: Server-side map generation from postGIS database

texts:
I have loaded some OSM data into a postGIS database, using imposm.
Now using this data base, I computed a value associated to every point of the map (could be anything, real-estate average price or yearly rainfall...).

Now I wish to draw the map and to represent that information on it with colors. Ideally, I need to generate both a neat classical jpg image and also tiles to be displayed on a website. I'm looking for a nice quality output.

I'm not looking for a GUI tool, I would like it to run on a server.

I had a look at Mapnik, I saw GeoServer and QGIS mentioned as well, but being new to the GIS subject I feel a bit confused, as I thought coloring a map would be a very vanilla task, but it seems not.

What tool would you please recommend ?


texts: Fixing "Error: creation of data source failed (OGR error:Failed to create directory .shp for shapefile datastore.)"

texts:
I am trying to convert from decimal degrees to meters in QGIS. I have a shapefile that is currently in WGS84/UTM zone 19S. I have tried the save as - browse - WGS84 EPSG:4326 but this appears:

Export to vector file failed. Error: creation of data source failed
(OGR error:Failed to create  directory .shp for shapefile datastore.)

Does anybody have any ideas? Or any other way of doing it?
I am using QGIS 2.2.0


texts: QGIS - Determine percentage of coverage of buffer for census tracts

texts: I have a census tract shapefile for the state of South Dakota. I have providers plotted with a 40 mile radius buffer. I am trying to determine what percent of each census tract is covered by any part of the buffer. How can I do this?


texts: I'm trying to locate or find out if there is a standard style for municipal infrastructure.  Specifically Canadian (Provincial BC, AB) standards for symbolizing water, sewer, valves, hydrants etc.


texts: Is geotools transformation of EPSG:3785 (Web mercator) to EPSG:4326 wrong?

texts:
Using postgres (9.3.6) and postgis:

 select st_astext(st_transform(ST_GeomFromEWKT('SRID=3785;POINT( 7550310.77 4139679.02 )'), 4326));

This returns POINT(67.8255956454323 34.8217133646348)

Using scala and geotools (11.2 but I tried 13.0 as well):

val c = new Coordinate(7550310.77, 4139679.02)
val hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, true)
val crsAuthFactory = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", hints)
val wgs84 = crsAuthFactory.createCoordinateReferenceSystem("EPSG:4326")
val webCRS = crsAuthFactory.createCoordinateReferenceSystem("EPSG:3785")
val webToWorldTransform = CRS.findMathTransform(webCRS, wgs84)
JTS.transform( c, null, webToWorldTransform )

This returns (67.82559564543232, 35.00213150097441)

The latitude seems to be shifted using geotools. Any advice? Am I missing something?


texts:
I have kml files that I can see on google eart and street view, the problem is that when I load the same data on the google earth app, I see them in 3d and 2d but not when I go down in street view while before I saw them. Where am I wrong?
enter image description here


texts: Using multiple workspace environments ArcGIS Arcpy

texts:
I have 1 folder with shapefiles and a geodatebase. I'm trying to export the shapefiles and geodatebase feature datasets but I need to reference 2 Workspace Environments in order for the script to read both destinations.

Error Message:
Runtime error  Traceback (most recent call last):   File "", line 68, in    File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\conversion.py", line 1811, in FeatureClassToGeodatabase     raise e ExecuteError: ERROR 000732: Input Features: Dataset 'ODORANT TEST POINTS.shp';COMPONENTS.shp;FARM_TAPS.shp;LINE_MARKERS.shp;MAIN_LINES.shp;METERS.shp;SERVICE_LINES.shp;STATIONS.shp does not exist or is not supported

  # Read (2) Workspace environments
  67. for ws in workspaces:
  68.     arcpy.env.workspace = ws

Code:

# Import modules
import arcpy

arcpy.env.overwriteOutput = True

# Set environment settings
workspaces = [
    'M:/EIS_DRAFTS/EIS/Company SHP Files'
    'M:/EIS_DRAFTS/EIS/Company SHP Files/Company UPDM.gdb'
]

# Set local variables
in_dataset_P_Integrity = ['P_Anomaly', 'P_AnomalyGroup',
'P_ConsequenceSegment',
'P_ControlPoint',
'P_CouldAffectSegment',
'P_DASurveyReadings',
'P_DocumentPoint',
'P_DOTClass',
'P_Elevation',
'P_GasOdorReport',
'P_HighConsequenceArea',
'P_ILIGroundRefMarkers',
'P_ILIInspectionRange',
'P_ILISurveyGroup',
'P_ILISurveyReadings',
'P_InlineInspection',
'P_InspectionNote',
'P_InspectionRange',
'P_LeakSurvey_AreaSurvey',
'P_LeakSurvey_Atmosphere',
'P_LeakSurvey_BareSteel',
'P_LeakSurvey_Business',
'P_LeakSurvey_CastIron',
'P_MAOPCalcRange',
'P_OperatingPressureRange',
'P_PipeCrossing',
'P_PipeExposure',
'P_ServiceInterruptionEvent',
'P_StationSeries',
'P_TestPressureRange']

in_dataset_P_PipeSystem = ['P_CompressorStation', 'P_ControllableFitting',
'P_CPAnode',
'P_CPBondJunction',
'P_CPBondWire',
'P_CPRectifier',
'P_CPRectifierCable',
'P_CPTestPoint',
'P_DehydrationEquip',
'P_Drip',
'P_ExcessFlowValve',
'P_GasLamp',
'P_GatheringFieldPipe',
'P_LineHeater',
'P_MeterSetting',
'P_NonControllableFitting',
'P_Odorizer',
'P_PigStructure',
'P_Pipes',
'P_PipeSystem_Net_Junctions',
'P_PressureMonitoringDevice',
'P_PumpStation',
'P_RegulatorStation',
'P_ReliefValve',
'P_RiverCrossing',
'P_RuralTap',
'P_Scrubber',
'P_Service',
'P_Strainer',
'P_Tank',
'P_TownBorderStation',
'P_Valve',
'P_Well']

in_dataset_P_EverythingElse = ['ODORANT TEST POINTS.shp',
'COMPONENTS.shp',
'FARM_TAPS.shp',
'LINE_MARKERS.shp',
'MAIN_LINES.shp',
'METERS.shp',
'SERVICE_LINES.shp',
'STATIONS.shp']

out_dataset_P_Integrity = 'C:/EnterpriseFolder/Company.gdb/P_Integrity'
out_dataset_P_PipeSystem = 'C:/EnterpriseFolder/Company.gdb/P_PipeSystem'
out_dataset_P_EverythingElse = 'C:/EnterpriseFolder/Company.gdb'

try:
    for ws in workspaces:
        arcpy.env.workspace = ws
        # Use ListFeatureClasses to generate a list of inputs
        for infc in arcpy.ListFeatureClasses():

        arcpy.env.outputMFlag = "Disabled"
        arcpy.env.outputZFlag = "Disabled"

        #Execute FeatureClassToGeodatabase
        arcpy.FeatureClassToGeodatabase_conversion(in_dataset_P_Integrity,out_dataset_P_Integrity)


texts: Where is information on What's Coming in ArcGIS 10.3 for Desktop?

texts:
Every year I need to make a case to my employer as to why it is essential we continue to pay maintenance.

To justify this investment, I need to have some idea about what's coming so I began googling what's new in ArcGIS 10.3 Desktop and I can't find any concrete information.

Does anyone know where I can find this info?

I understand ArcGIS Pro comes bundled with 10.3, but at this time I am only interested in new features, fixes and advancements in ArcGIS 10.3 for Desktop and its extensions.


texts:
I have updated PostGIS to 2.4.3. Using with Postgresql 9.6 on Ubuntu 14. When I try to select ST_AsMvt(example_geometry) from example_table, I get:

 ERROR:  Missing libprotobuf-c
 ********** Error **********

 ERROR: Missing libprotobuf-c
 SQL state: XX000

As it it mentioned on PostGIS Installation:

  To enable ST_AsMVT protobuf-c library (for usage) and the protoc-c
  compiler (for building) are required. Also, pkg-config is required to
  verify the correct minimum version of protobuf-c.

If I execute protoc --version, I get libprotoc 3.5.0.

libprotobuf-c0, libprotobuf-c0-dev, libprotobuf-dev are already installed.

Tried sudo apt-get install --reinstall postgis but didn't help.

Related issues:

ST_AsMVT Centos 7 Missing libprotobuf-c
Postgis ST_AsMVT() return missing libprotobuf-c on MacOSX High Sierra

How to solve this problem?


texts: Exporting BNG coordinates from MapInfo (27700) without rounding

texts:
I have a TAB file in ESPG 27700 projection and MapInfo is rounding Y coordinate values upon export. It also occurs when viewing the table (although selecting the value shows its 'true' form)

Upon export to .csv the value stays at 53.5834. I have tried creating a character field and updating it with the MIDY value and this doesn't work either.

I am currently experimenting with generating a general text direction assigned to a polyline by first splitting it to lines and using the difference between the X and Y coordinates of the centreline of each segment and 'counting' the predominant X and Y directions (disregarding values <=5m, the 'NEG' values) with the aim of possible re-creating my steps in MapBasic:


texts: PostGis Scaling for Polygons at a fixed center location

texts:
I need to enlarge/expand a 2D polygon at its fixed location (lat,lon) with a constant scale factor in order to retrieve the surrounding area with the same shape as the origin polygon.

What PostGIS function can I use and is this even possible to achieve?

Principally, I'm looking for a function similar to PostGISs ST_Scale(), which only works for geometry data, which has no center.

Can you assist with the necessary transformations that come with this issue?

My Idea to solve the problem. I'm sure, that it can be formatted better:

SELECT
st_astext(st_centroid(way)) as polygoncenter,
ST_AsText(way) as originalpolygon,
ST_AsText(
    ST_Translate(way, -ST_X(ST_AsText(ST_Centroid(way))), -ST_Y(ST_AsText(ST_Centroid(way))))
    ) as translatedtoorigin,
ST_AsText(
    ST_SCALE(ST_Translate(way, -ST_X(ST_AsText(ST_Centroid(way))), -ST_Y(ST_AsText(ST_Centroid(way)))),2,2)
    ) as scaledatorigin,
ST_AsText(
    ST_Translate(ST_SCALE(ST_Translate(way, -ST_X(ST_AsText(ST_Centroid(way))), -ST_Y(ST_AsText(ST_Centroid(way)))),2,2),ST_X(ST_AsText(ST_Centroid(way))), ST_Y(ST_AsText(ST_Centroid(way))))
    ) as scaledatoriginalposition
FROM de_polygon LIMIT 1000

The last parameter scaledatoriginalposition is the value I want to retrieve.
Isn't there any shorter way or given functionality by Postgis to solve this problem?


texts: Conditional Combine Rasters and Track which Raster pixel used

texts:
I have 6 rasters with pixel values ranging from 0 to 1. I want to extract the highest pixel values from each raster to create a new raster, but also have a way to show which raster the highest value came from.

For example, pixel N has a value of 0.1 in Rast1, 0.5 in Rast2, 0.4 in Rast3.
I would like to keep the value from Rast2 in the output raster and somehow keep track that this value came from Rast2 not Rast 1 or 3. Ultimately I would like to symbolize the raster by which original file it came from.

Is this possible?


texts: Reduce WFS response size by reducing feature detail

texts: Using GeoServer 2.5.4 is there a way to request lower detail features? I currently have a 4.75MB response and each coordinate pair is extremely precise. Since the .shp providing the features may change occasionally I don't want to spend time optimising the file itself. I'd rather configure GeoServer or my WFS request to only return the coordinates to within 1m precision, or to only return those necessary according to some optimisation algorithm such as Douglas-Peucker. Is that possible?


texts:
I downloaded a Nigeria shapefile and loaded it via the add vector layer where I picked out the three .shp files within the zip file and displayed them on the layers panel. I can see the shapes displayed on the canvas.

On the canvas there are no state names or local government names so I thought those attributes may be in a .csv file, included in the zip file. I tried to add it as a delimited text file but the file doesn't contain any coordinates.

Nothing is displayed on the canvas, no state names or anything. How do I make the attributes within the .csv file to be displayed on the canvas?

I'd like to mention that I'm a beginner in GIS and I just downloaded the QGIS software only a couple of hours ago.


texts: pgr_ksp() returns other columns from edge table

texts:
I'm willing to optimise my queries to make it much faster.

In my Edge table, I have different information for each edges (distance, speed, duration, ...).
The problem is so far my pgr_ksp() function only can return these data: seq, path_id, path_seq, node, edge, cost and agg_cost.

So far what I do to get the complete data needed, I need to add subqueries to return the complete data depending on the current edge of the sequence such as bellow:

SELECT seq, path_id AS route, path_seq, node, edge,
       cost as duration,
       (SELECT distance FROM edges WHERE id=edge) as distance,
       (SELECT speed FROM edges WHERE id=edge) as speed,
       agg_cost
FROM pgr_ksp('SELECT id,
              source::integer,
              target::integer,
              duration::double precision AS cost,
              reverse_cost::double precision
       FROM edges
       ORDER BY id',
       source_id, target_id, alternatives, true, false);

How could I do to avoid subqueries and therefore improve query speed?


texts: QGIS: Georeferencing a scanned map with known LCC projection and lat/long grid

texts:
I looking for advice or recommanded methods to georeference scanned maps for which I know the projection and some coordinates with lat/long dat in degrees

In my test case, projection is:

lcc: Lambert Conformal Conic  - 93 / IGNF:LAMB93)

LAMB93 requires coordinates in meters.

The map contains several cross-hairs marking and I have marked some places with GPS coordinates.

I have started to georeference the raster(scanned map) with QGIS Plugin using these known points with lat/long coordinates in degrees on the map. I used the tutorials here: Georeferencing Topo Sheets and Scanned Map and here: 15.2. Lesson: Georeferencing a Map )

My concern is that when importing the map in QGIS and superposing on other maps, I have an offset (and/or due to a slight rotation).

I assume that by giving coordinates in degrees, QGIS imports the whole map as "a WGS84 raster". I did not find a way to specify target SRC with a lcc projection.

I tried to used RGF93 (EPSG 4171), but it ended also with an offset.

My aim is to find a method to import the scanned maps within the original projection (SCR) and correlating the control points with
as they have been done

NB: I managed to reduced the offset by using transformation distortion of the imported map and multiplying the correlation points, but I am not satisfied with the result and extra work required.

Notes:

QGIS  version 3.4 used

-LAMB93 proj file:

Proj4: +proj=lcc +nadgrids=@null +a=6378137.0000 +rf=298.2572221010000 +lat_0=46.500000000 +lon_0=3.000000000 +lat_1=44.000000000 +lat_2=49.000000000 +x_0=700000.000 +y_0=6600000.000 +units=m +no_defs

-RFG93 proj file:

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

QGis Georeferencer, Transformation Window Plugin Window


texts: Weighted averages: accounting for grid cell surface variation in proportion to the cosine of its latitude

texts:
Non-geographer learning GIS here. I'm using ArcGIS Desktop v. 10.2.2.

I have a 30-arc second Terrain Ruggedness Index (TRI) raster file, and I want to measure the average ruggedness within each polygon of shapefile I have (not public). The TRI data is available here.

The authors of these data point out that "it is important to take into account that the sea-level surface that corresponds to a 30 by 30 arcsecond cell varies in proportion to the cosine of its latitude." They provide a data file measuring the surface area of each cell, in square metres, which must be used to weight the TRI when averaging across polygons.

Questions:

Instead, can I re-project the raster data to an equal area projection and then calculate an unweighted average raster value? In other words, does an equal area projection account for this surface variation in proportion to the cosine of its latitude?
If the weighted average raster value route is a must, is the appropriate method to use the Times tool with the TRI raster and cell area raster, then calculate the average raster value of each polygon with Zonal Statistics?


texts: GDAL - Register Pixel Function for Command-Line Utilities

texts:
I'm new to GDAL and am just learning the inner-workings of C++ compilation at the moment so please bear with me.
I am trying to define a pixel function to produce a VRT from input TIFs, with a band being some computation of other bands, which I can translate into a TIF. Let's say something simple like the min of each band.
The top voted answer here suggests that these functions cannot exist outside the program where they are registered. There is a nice post from the GDAL mailing list which summarises the process which suggests otherwise however. This post indicates that I should be able to build a shared library, and define the GDAL_DRIVER_PATH environment variable to point at the directory containing this. This looks promising, since now a call to gdal_translate on the command-line shows that it's trying to access this library but can't find the GDALRegisterMe or GDALRegister_pixelFunctions function (my library is called gdal_pixelFunctions.so)
$ gdal_translate testVrt.vrt testVrt.tif
ERROR 1: build/gdal_pixelFunctions.so: undefined symbol: GDALRegisterMe
ERROR 1: build/gdal_pixelFunctions.so: undefined symbol: GDALRegister_pixelFunctions
Input file size is 1000, 1000
0ERROR 5: VRTDerivedRasterBand::IRasterIO:Derived band pixel function 'OneFunction' not registered.

My src looks like
CPLErr OneFunction(
    void **papoSources,
    int nSources,
    void *pData,
    int nXSize,
    int nYSize,
    GDALDataType eSrcType,
    GDALDataType eBufType,
    int nPixelSpace,
    int nLineSpace
)
{
    // trivial implementation
}

void GDALRegisterMe()
{
    GDALAddDerivedBandPixelFunc("OneFunction", OneFunction);
}

Attempting this on Arch Linux if that helps.
EDIT:
In case this is helpful - This link shows how GDAL tries to load drivers from shared libraries, searching first for the GDALRegister_pixelFunction() function (in this case) and then GDALRegisterMe() if not found, but my shared library seems to be failing to provide these for some reason.


texts: Cutting polygons in same isochrone map layer?

texts:
I've created some travel time isochrones using the 'OSM Tools' plugin in QGIS. I now need to remove overlapping sections of the isochrone layer so that I'm left with discrete polygons for each isochrone (in this case travel time to a central point). This is so I can calculate the proportion of a 'district' that is within x and y and z, mins of the city centre. So for example here are my isochrones

If I now remove the 'shortest travel time' (yellow polygon) it's clear that all other polygons depicting travel times in excess of the yellow polygon are beneath:

I would like to cookie cut out the longer travel time polygons so that none of the individual isochrones overlap i.e. my layer looks like this (with the yellow isochrone hidden):

I had thought that the difference or symmetrical differnece tool would do the trick but no joy. Is there an easy way to cut/clip polygons on the same layer from within QGIS?


texts: How to fix "Detected out of bounds geographic extent"

texts:
I hope someone can help me to solve this problem.
i want to make a map with TileMill v.01, when i try to add layer using shapefile (i use QGIS to edit *.shp file) or postGIS i always get error message :  Detected out of bounds geographic extent (289148.0964570158,9020963.6199987,180,85.051) for layer 'badung'. Please ensure that the SRS for this layer is correct. Its native extent is '289148.0964570158,9020963.6199987,307432.8551368474,9088475.5022231'

i think if i change the geom property in postGIS using this WGS84 (4326) will fix this problem, apparently not. could someone give me the way to fix this problem ?


texts: How to extract list of XY coordinates found within a polygon?

texts:
I have multiple polygon shapes covering the UK. How can I create an Excel spreadsheet of all 6 (sometimes 7) figure XY British National Grid northings/eastings that fall within these polygons?

I don't want just the XYs of the polygon verticies, but all coordinates intersecting the polygon.


texts: How to use an ArcGIS Server cache directory with multiple (differently named) map services?

texts:
We have a map service named 'X' that created and uses tiles in a cache directory 'Folder_X'.

I'd like to create a second map service named 'Y' that uses the same cache directory - that is, it serves up the same tiles by a different service name.

So far, I see that when creating a service in ArcCatalog, it insists on using a cache directory named to match the service name, which of course isn't what I want. I see that there is a 'Folder_X'.cfg text file in the ...ArcGIS\server\user\cfg\Folder directory that I can edit to point to a different folder, but it seems to keep getting overwritten with the default value.

Ideas?


texts:
My question is about reverse geocoding addresses over a large area from a set of points or polygons in ArcMap.

I know from previous experience that Alteryx with the Spatial Match tool allows one to import Google Earth KML files, after which Spatial Match will return addresses. The KML can contain either a polygon (which returns addresses within the polygon), or a set of points (which return addresses at those points).

In ArcMap, I can use the Address Inspector tool to see the data I want, but can't find a way to batch the data out over an area. Setting the search distance doesn't seem to help either. I currently have ESRI World Geocoder and MGRS address locators. Also, it seems I can only import the Google Earth KMLs as a layer file.

Ultimately, how do I do reverse geocode from Google Earth KML on a large scale?


texts:
I tried to run gdallocationinfo from my FWTools shell and got a message saying: "gdallocationinfo is not recognized as an internal or external command, operable program or batch file". So I checked the installation folder and didn't see any file with this name.

So I am guessing that something could be missing.
I've also tried to search for the file on Google but without success.
I am using FWTools 2.4.7 installed on Windows 7.

Any idea on how to fix this?

Thanks


texts: osm data import failed with osm2pgsql and windows

texts:
I am trying to import a .pbf file into my postgres/postgis database. I am using the following command line for importing the osm file

  C:\Users\Toton>osm2pgsql -c -d osmgb -U postgres -H localhost -S
  C:\osm2pgsql\os m2pgsql-master\default.style
  E:\osmgb\great-britain-latest.osm.pbf

Its returning the following error

  release notes: 'Windows version built by Dominik Perpeet
  (http://www.customdebug.com/osm2pgsql/index.html)'  Using projection
  SRS 900913 (Spherical Mercator)  Error reading style file line 151
  (fields=4) flag 'phstore' is invalid in non-hstore mode Error
  occurred, cleaning up


texts: Displaying road outlines but not segments i.e. no lines at intersections?

texts:
I have the roads that I want to just display at outlines but not segments (I.e not disconnected from each road).

I am trying dissolve but its not working. Attached is a link to a package file of the data. Please help!  https://www.dropbox.com/s/55hqr9471v9qxwf/Road%20Test%20package.mpk


texts: How can I build CRS string in Trasverse Mercator to increase distance measurement accuracy?

texts:
I've searched a lot in the forum finding several answers to part of my problem. In brief I need to make some computation keeping in account the distance from every point along different angles until a distance of 25 km. The resolution of raster from which I'll extract data change as step function of distance from the point analyzed. It goes from 1m in a narrow zone (1 km) to around 90m (until 25 km). The zone is Europe. This give an indication on the accuracy I need.
Initially I've been starting using UTM implementing in the code the simple capability to select the correct UTM zone. To minimize distortion when I get closer to the UTM zone limit (or higher lat) I'd like to change approach implementing a TM projection in Elliptical form centered in the area I'm analyzing that should means define as central meridian of grid origin the mean values of LON  of the area considered (solution partially suggested by this answer of @whuber in Measuring distance in spherical Mercator vs zoned UTM).

How can I build the projection string in a parametric way?

My main problems, cause my limited knowledge in this topic, is about some parameters as +k +x_0 +y_0 (Scale factor at natural origin, False easting and False northing).
The CRS string should appear in this way if I'm not wrong:

    "+proj=tmerc +datum=WGS84 +lat_0=0 +lon_0=b +k=? +x_0=? +y_0=? +units=m
     +no_defs +ellps=WGS84 +towgs84=0,0,0"

with b the mean lon of the area analyzed.
For my knowledge false easting and northing should only affect the sign of coordinates but +k?
Can I set +k as 1?

Is the approach correct? It is better to choice a different map projection in my case?


texts: How close can a commonly used GPS device come within actual location of target

texts: With resepct to meters/feet how close will a GPS device get me, say if device (like iPad), is "common". Say, could I upload a particular coordinate that has two-tree houses in my back yard and discern one from the other? I want to build an app which discerns between plots separated by just a few feet.  Possible and plausible?


texts:
I am interacting with a google map using openlayers3. On pinch zoom, the map throws an exception : Uncaught TypeError: Cannot read property 'x' of undefined

I am using Google Maps API 3.2 and OpenLayers 3.0.0.

EDIT

I notice the map is losing zoom property on pinch zoom:
googleMap.getZoom() now returns undefined. Thus the exception is thrown.


texts: How to set up a threshold for Delaunay triangulation?

texts:
I've been working Delaunay triangulation on a set of random points inside a group of mostly linear areas (in fact, buffers around rivers).

The problem is that while the triangulation works correctly inside these areas, it also produce a lot of triangles between the areas. Expectedly.

Now, I need to delete all this extra data manually.

Does anyone think it would be possible to define a threshold, a maximum distance, for the Delaunay algorithm?

That way, it would only produce triangles inside each area.

Any idea how I could do that? I'm using QGIS or GRASS.


texts: Copying specific folder into subfolders with ArcPy?

texts:
I'm working on Arcgis 10.5 and I am trying to achieve this:

On my script, I want to make folders and subfolder, and also I want to copy one specific folder into several subfolders at the same time. My script is working well to make the folders and subfolder.

However, I get  the next error when the script tried to copy the folder I want :

  TypeError: cannot concatenate 'str' and 'Result' objects.

This my script:

I can't see where is the mistake so can somebody can give me a hand?

import arcpy
import os
import shutil
#Workspaces
arcpy.env.workspace = r'F:\FOLDER'
Folder_to_copy = r'F:\FOLDER\FOLDER1\FOLDER2\BASE_CHIH_MZA_UTM13_2018' #<--I want to copy this folder into the subfolders Sub1,Sub2,Sub3,Sub4

#Folder to make
Sub1_Folders = ['Sub1FolderA','Sub1FolderB','Sub1FolderC']
Sub2_Folders = ['Sub1','Sub2','Sub3','Sub4'] #<--The Especific folders where I want to copy my folder

#Result
Out_putfolder = r'F:\FOLDER\OutPutFolder'

#Field dbf to make subfolder
Tabla_dbf = r'F:\FOLDER\BRIGADISTA.dbf'
Field1 = 'BRIGADISTA'

Brigadistas = [row[0] for row in arcpy.da.SearchCursor(Tabla_dbf,Field1)]
Unicos_Briga = set(Brigadistas)

#Make folder and subfolders
for f1  in Sub1_Folders:
    Sub1folder = arcpy.CreateFolder_management(Out_putfolder,f1)
    for f2 in Sub2_Folders:
        Sub2folders = arcpy.CreateFolder_management(Sub1folder,os.path.join(f2))
        for tb in Unicos_Briga:
            tabladbf = arcpy.CreateFolder_management(Sub2folders,os.path.join(tb))
            Copyfolder = shutil.copytree(Folder_to_copy,os.path.join(tabladbf))

print 'Done'

This is the error show up

Runtime error
Traceback (most recent call last):
  File "<string>", line 29, in <module>
  File "C:\Python27\ArcGIS10.5\Lib\ntpath.py", line 90, in join
    return result_drive + result_path
TypeError: cannot concatenate 'str' and 'Result' objects


texts: Trouble clipping RGB raster in QGIS: Layer not created

texts:
I am brand new to QGIS, and am having trouble clipping a raster to a shape file.
Starting with the three-band red-green-blue (RGB), named 'ortho':

...I created a shape file to extract the plots that I wanted to, named 'field_boundary':

I am using the GDAL plugin in order to do this. I went to the top taskbar and chose 'Raster' -> 'Extraction' -> 'Clip Raster by Mask Layer.' The following parameters were used to run the command:

I am trying to generate a new raster, as the Clip tool would do if using ArcGIS. However, I get the following error:

What can I do to get this to clip to a new raster?


texts: converting UTM - Decimal Degrees with javaScript or using a web service

texts:
I am creating a wildlife survey form that has a field where users enter location. Users will be able to enter location using a map (i.e. draggable pushpin) or by manually entering lat/lng. My clients are also requesting an option for entering UTM coordinates.

If a user does enter UTM coords I would like to convert these before pushing to the database. I've seen a couple websites that have utilities for converting UTM to decimal degrees.

Does anyone know of a javaScript lib that will make these conversion given zone and easting-northing? I also saw this thread

Convert geographic coordinates to UTM zone 29N with Javascript

Does ESRI have a web service that I could pass UTM coords to and get lat/lng back?


texts:
I have some data in a database and would like to show them on Leaflet, but unfortunately I cannot.

My data look like:

 Columnname:Val1
 POINT (-36.389819 174.124343)

and my code:

var map = L.map('map').setView([-36.799907, 174.708520], 11);
L.tileLayer('http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}', {
         attribution: 'Tiles &copy; Esri &mdash; Source: Esri, DeLorme, NAVTEQ, USGS, Intermap, iPC, NRCAN, Esri Japan, METI, Esri China (Hong Kong), Esri (Thailand), TomTom, 2012'
    }).addTo(map);

var mywms = L.Point.wms("http://localhost:9090/geoserver/localhost/wms", {
                 layers: 'localhost:DP_View',
                format: 'image/png',
                 transparent: true,
               version: '1.1.0'
    });
mywms.addTo(map);

Seems L.Point does not have a method for WMS, not sure how I can show these data on the map with a popup. I am able to see this data as PreviewLayer in GeoServer.

What I have tried so far:

 Changed WMS style as suggested

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld
StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
 <!-- a Named Layer is the basic building block of an SLD document -->
<NamedLayer>
<Name>default_point</Name>
<UserStyle>
<!-- Styles can have names, titles and abstracts -->
  <Title>Default Point</Title>
  <Abstract>A sample style that draws a point</Abstract>
  <!-- FeatureTypeStyles describe how to render different features -->
  <!-- A FeatureTypeStyle for rendering points -->
  <FeatureTypeStyle>
 <Rule>
   <PointSymbolizer>
     <Graphic>
       <Mark>
         <WellKnownName>square</WellKnownName>
         <Fill>
           <CssParameter name="fill">#009900</CssParameter>
         </Fill>
       </Mark>
       <Size>12</Size>
       <Rotation>45</Rotation>
     </Graphic>
   </PointSymbolizer>
   </Rule>
  </FeatureTypeStyle>
</UserStyle>
 </NamedLayer>
</StyledLayerDescriptor>

The output

And the Spatial Results from MSSQL


texts: How can I delete multiple field attribute values from an attribute table?

texts:
I have a Centerlines feature class to which I transferred attributes from an AddressPoints feature class, and now each centerline has the correct street name but also a street number which I don't need since it is a polyline.
I'd like to clear those values from the table, without deleting the street name field.


texts: Bringing attributes of multiple polygons into other polygons?

texts: I need to absorb multible attributes included in block areas into a number of ploygons. In other words, I need to bring block level data into a number of overlaying polygons. The polygons will absord all the data in the smaller polygons. They can lose their data so long as their spatial location remain the same. Below is an image of the project. The blue zones will be absorbed into the Orange zones.


texts:
My goal is to test and compare the performance of 3D tiles VS i3S in :

storing large volumes of data (volumetry).
displaying (data display time)

I have currently 2 larges 3D tiles files. For information, I work on Windows OS.
Is there a way to convert those two 3D tiles to i3S format ?


texts: C++ GDAL API: Add more than one rasters to VRT

texts:
I want to add more than one tiff file into a vrt file, Well Based on GDAL documentations it is possible to create a vrt using this code

GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
GDALDataset *poSrcDS, *poVRTDS;
poSrcDS = (GDALDataset *) GDALOpenShared( "utm.tif", GA_ReadOnly );
poVRTDS = poDriver->CreateCopy( "wrk.vrt", poSrcDS, FALSE, NULL, NULL, NULL );
GDALClose((GDALDatasetH) poVRTDS);
GDALClose((GDALDatasetH) poSrcDS);

but it only adds utm.tif. what If I want to add another set of images into this vrt file? can I use poDriver->CreateCopy again and it appends the images into the source vrt?


texts:
I am in QGIS.
I have a raster layer of mountains and a lake.
I have a vector layer (shapefile) of the lake. It sits nicely where it supposed to be over the raster level.
I am going to convert this height map to an STL so that I can cnc this out of wood, so I'd like the actual lake to be something like 200 feet lower than it shows on the raster layer so that I can fill it with blue epoxy.
I have read a lot of clipping and everything else, and I am lost.


texts: Displaying bounding box polygon from CSV using ArcGIS Desktop?

texts:
I have a data file (.csv) with a minimum and maximum longitude and latitude. These 4 fields (roughly a few thousand rows) make up a bounding box. However, I am unsure of how to get this data to be displayed as a polygon. When I try to display data, I can only do points(X,Y). I need to create an Xmin, Xmax, Ymin, Ymax for each row of data.

I have brief experience with ArcPy, but would like to avoid unless that is your only solution. Also, I have this file also as a .JSON. I tried using the convert .JSON to features tool, but there was an error saying unable to parse when run.


texts: Exporting Adobe Illustrator (AI) file layer from QGIS?

texts:
How do I export Adobe Illustrator (AI) layers from QGIS?

I have tried all the "save as" files from layer and non of them worked.


texts: Opening shapefiles without prj-file (one of Gauss–Krüger CS) in QGIS

texts:
I'm trying to open a shapefile in Gauss–Krüger coordinate system, which comes without any .prj-files, in QGIS. These files are not really up to date but they're the only data i can get.
Is it possible to project them on any world map in the right position by using the "OpenLayersPlugin"?
The other part of the data is more updatet so its CRS is ETRS 89 / UTM 32N. Is there a way to project both on the same map for the same time?
When I insert the Gauss–Krüger data into QGIS, it says that the CRS of the Gauss–Krüger data is undefined and preset to EPSG:4326. What does that mean?
So I need a way to show both, old and new data, beside each other on the right position in the map (Google Satellite from OpenLayersPlugin).


texts:
I have a geotiff that I have been working with in tilemill and I would like to select portions (buildings) in the map and create new layers as geojson objects. My issue is I don't know how to get the coordinates to create the objects!

http://chimera.labs.oreilly.com/books/1230000000345/ch12.html#_choropleth

this is a link to what I'm trying to do with states I just need the step where the coordinates are gotten because it doesn't seem to be in this document, they just give a shout to the guy who defined them. If anyone has done something similar to this, please give me a shout!

TO CLARIFY

I am not looking for the coordinates for the 50 states, I am trying to create geojson objects for a FICTIONAL map.


texts: Load the raster image which is in TIF format using QGIS API

texts:
I am quite new inthe field of GIS.
I am interested in is to load the raster image which is in TIF format using QGIS API

I have followed the steps given the in the follwing link to load shape file
enter link description here
http://geospatialdesktop.com/2009/02/creating_a_standalone_gis_application_1/ and it works fine

To load the rester image i have followed the steps given int he link
http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html

I have python 2.7.6 ,GDAL 1.10.1 installed with qgis2.0 installed in my system

code snippet

import osgeo.gdal
from  gdalconst import *
from PyQt4.QtCore import QFileInfo,QSettings
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem

file = QFileDialog.getOpenFileName(self,
    "Open TIF Image", ".", "Tagged image (*.tif)")
fileInfo = QFileInfo(file)
fileName=fileInfo.fileName()
baseName = fileInfo.baseName()

layer = QgsRasterLayer(fileName, baseName)

but the image is not loaded and i got the following error in the terminal

ERROR 4: `mypic.tif' does not exist in the file system,
and is not recognised as a supported dataset name.

Layer failed to load!
QObject::connect: Cannot connect (null)::repaintRequested() to QgsMapCanvas::refresh()
QObject::connect: Cannot connect (null)::screenUpdateRequested() to QgsMapCanvas::updateMap()


texts: Concerning "Error when using r.series in a python script"

texts:
Based on the "Error when using r.series in a python script" discussion (Error when using r.series in a python script) and based upon the solution given by Joseph I tried to make the same but by using a loades set of rasters and saving the result as temporary file. I made the following script but it doesn't work to me. Any Idea to fix it?

##RasterMaxFinder=name
##1=multiple raster
##maximum=output raster

import glob, os, processing
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsRasterLayer, QgsRectangle

raster = processing.getObject(1)

list = []
extent = QgsRectangle()
extent.setMinimal()

for raster in 1:
    fileInfo = QFileInfo(raster)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(raster, baseName)

    # Combine raster layers to list
    list.append(rlayer)

    # Combine raster extents
    extent.combineExtentWith(rlayer.extent())

# Get extent
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

# Run algorithm
processing.runalg('grass:r.series', list ,False,6,'-10000000000,10000000000',"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),0.0, maximum)


texts: Permanently replace/update server/service in each ArcSDE Layer in MXD?

texts:
Using the mxd.findAndReplaceWorkspacePaths method, if and MXD is updated and the path the the .sde connection file moves again all of the data sources will break again.

Is there a method to permanently update each SDE layer in an MXD Table of Contents to permanently store the SDE Server and Service? Ie. Just replace the SERVER and INSTANCE for each SDE layer in an MXD using python or ArcObjects in Pyhton?


texts: Publishing subtypes on ArcGIS Online for Collector?

texts:
I am kind of new to ArcGIS online and I just publishes a feature hosted map service after I have fixed my data with subtypes and subtypes-domains. After publishing to agol I discovered that all my domains are not working.

I need this for collector and I had intentionally created the domains based on subtypes in order to restrict the data entry based on the subtype.


texts:
I'm using Arc 10.1. I am creating a parcel inventory where I would like to have a photo of the property attached as an attribute and view-able with an HTML popup. I have a polygon feature class built in my geodatabase. I created a Photo Field with data type as raster. I can load a photo while in an edit session of ArcMap but when I chose save edits I get an SQL error that looks something like this:

"An invalid SQL statement was used [fras_blk_MyPropertyTable]
 An invalid SQL statement was used [SELECT * fras_blk_MyPropertyTable WHERE rasterband_id in()]"

Any thoughts?


texts: OpenLayers or OGCServer request optimization

texts: I have an OpenLayers-OGCServer-Mapnik setup through a WMS layer. I noticed when i use the openlayers UI in a quick manner (fast zoom with mouse scroll, a couple of consecutive pans with mouse) the response from WMS slows down, because its obviously processing all the tiles between the previous "static" state and the last one. Is there a switch in openlayers or perhaps mapnik ogcserver to optimize this behavior by either postponing requests until the user pauses for a moment on the same location (clientside) or canceling previous requests from the same IP on the OGCserver (serverside)


texts:
I want to migrate to Geonetwork 3.8.1 since 3.6 then I have taken some measurements

BACKUP my DB PostgreSQL
I change these files in order to point my new DB
              C:\Program Files\Apache Software Foundation\Tomcat 9.0\webapps\geonetwork\WEB-INF\config-db\jdbc
                  C:\Program Files\Apache Software Foundation\Tomcat 9.0\webapps\geonetwork\WEB-INF\config-node\srv

And finaly I replace the file metadata_data  by the older one in data directory

After doing these operations  the metadata records have been recovered. But when I try to display the full view mode it indicates that the metadata records are not found.

for the first full view mode i had this error

2019-10-11 09:54:26,841 ERROR [jeeves] - Error occurred within a transaction
javax.xml.transform.TransformerConfigurationException: Failed to compile stylesheet. 1 error detected.
    at net.sf.saxon.PreparedStylesheet.prepare(PreparedStylesheet.java:176)
    at net.sf.saxon.TransformerFactoryImpl.newTemplates(TransformerFactoryImpl.java:139)
    at de.fzi.dbs.xml.transform.CachingTransformerFactory.newTransformer(CachingTransformerFactory.java:120)
    at de.fzi.dbs.xml.transform.CachingTransformerFactory.newTransformer(CachingTransformerFactory.java:79)
    at org.fao.geonet.utils.Xml.transform(Xml.java:490)
    at org.fao.geonet.utils.Xml.transform(Xml.java:376)
    at org.fao.geonet.api.records.formatters.XsltFormatter.format(XsltFormatter.java:136)
    at org.fao.geonet.api.records.formatters.FormatterApi$FormatMetadata.call(FormatterApi.java:819)
    at org.fao.geonet.api.records.formatters.FormatterApi.getRecordFormattedBy(FormatterApi.java:330)
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)
    at java.lang.reflect.Method.invoke(Unknown Source)
    at org.springframework.web.method.support.InvocableHandlerMethod.doInvoke(InvocableHandlerMethod.java:221)
    at org.springframework.web.method.support.InvocableHandlerMethod.invokeForRequest(InvocableHandlerMethod.java:136)
    at ...

and for the second one is that :

2019-10-11 09:54:26,853 ERROR [jeeves] - Error occurred within a transaction
javax.xml.transform.TransformerConfigurationException: Failed to compile stylesheet. 1 error detected.
    at net.sf.saxon.PreparedStylesheet.prepare(PreparedStylesheet.java:176)
    at net.sf.saxon.TransformerFactoryImpl.newTemplates(TransformerFactoryImpl.java:139)
    at de.fzi.dbs.xml.transform.CachingTransformerFactory.newTransformer(CachingTransformerFactory.java:120)
    at de.fzi.dbs.xml.transform.CachingTransformerFactory.newTransformer(CachingTransformerFactory.java:79)
    at org.fao.geonet.utils.Xml.transform(Xml.java:490)
    at org.fao.geonet.utils.Xml.transform(Xml.java:376)
    at org.fao.geonet.api.records.formatters.XsltFormatter.format(XsltFormatter.java:136)
    at org.fao.geonet.api.records.formatters.FormatterApi$FormatMetadata.call(FormatterApi.java:819)
    at org.fao.geonet.api.records.formatters.FormatterApi.getRecordFormattedBy(FormatterApi.java:330)
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)
    at java.lang.reflect.Method.invoke(Unknown Source)
    at org.springframework.web.method.support.InvocableHandlerMethod.doInvoke(InvocableHandlerMethod.java:221)
    at org.springframework.web.method.support.InvocableHandlerMethod.invokeForRequest(InvocableHandlerMethod.java:136)
    at ...


texts: Using R, how do I rotate a GeoTiff file with multiple bands and save it?

texts:
My existing code (see below) only reads and rotates one band but I want all 360 of my bands to rotate.
library(raster)

your_raster <- raster("path/to/raster.tif")

rotated_raster <- rotate(your_raster)


texts: Reading Large Shapefile in R - alternatives to readOGR?

texts:
I am trying to read in a large shapefile using readOGR. The shapefile is
68837 features and 14 fields. It takes a considerably long time to read at >30min.  I understand this is a well documented issue with readOGR, but have not been able to find an adequate, updated solution.

How are folks dealing with reading large spatial datasets in R? Are there newer packages that are more efficient at reading large datasets?


texts: Saving routes from ArcGIS Network Analyst between set of points as shapefile in ArcGIS?

texts: Does anybody know how to  save the routes from ArcGIS Network Analyst between set of points as a shapefile in ArcGIS?


texts: Differences between mWm-2sr-1(cm-1)-1 and W/m^2 * sr * um

texts:
I'm confused about the use of this quantity. For example eumetsat

on pag. 16 radiance is equal to:

  Physical Units = Cal_Offset + (Cal_Slope x Level 1.5 Pixel Count)

expressed in mWm-2 sr-1 (cm-1)-1

and table in page 20 represents dynamic range in W/m^2*sr*um unit.

The question is: when I read seviri spectral channel from MSG2/3 and after the transformation above, do I have values expressed in first or second units?


texts: Mercator projection system grid problem. Grid false on QGIS but right on ArcGIS

texts:
I am actually trying to make a map similar to this, but for years 2019 and 2020.

The shapefile of the country I am using is in World Mercator (EPSG: 54004).
I tried on ArcGIS, but working on symbols transparency is really difficult. Here, longitude and latitude coordinates are well displayed on the grid (See image below).

So, for symbols, I wanted to use QGIS, but I have a problem actually. I would like to have that grid to tell how are longitude and latitude and to be able to precisely show where it is displayed. But on QGIS, latitude and longitude shown are pure nonsense (See image below).
.
I already tried to setup another projection system in grid parameters. And I tried to play a little bit with my projection parameters but nothing is working.
How can I repair that ?
I'll paste my "grid" setup and my map setup below.


texts: Finding given objects under both a polygon and a buffer

texts:
I'm a bit new to ArcGis 10.1, but know the basics.

I have two source layers, a polygon and a huge buffer. The polygon contains some of that buffer, but not all of it. There are numerous points inside that polygon (and outside), and I want to know how many of them are inside both the polygon AND the buffer.

Select by location only allows me to ask for one of the source layers, not both at the same time.  Any ideas of how to proceed ?


texts: Why do ECWs go missing when exporting MXD to PDF with ArcPy?

texts:
I'm trying to export some MXD maps to PDF with arcpy (ArcGIS 10.1). Everything works fine and dandy except for ECW images, which don't appear in the PDF. However, when I execute the same code from the ArcGIS python console, the ECW is exported. So if I execute this from an external script:

import arcpy
mxd = arcpy.mapping.MapDocument("path\\to\\test.mxd")
arcpy.mapping.ExportToPDF(mxd, "path\\to\\test.pdf")

I don't get the ECW layer in the PDF. But if I execute the same lines in the console, it works and the ECW layer is exported properly. Also, if I run

execfile('export_mxd_to_pdf.py')

from the ArcGIS console, I also get the ECW in the PDF. I don't have this problem with other image formats, like JPG.

What's going on here? What's different between the environments? How can I debug this?


texts: Creating dynamically multiple layers in ArcGIS for Server?

texts:
I am using ArcGIS Server 10.1 and Desktop 10.1 software.

I am publishing table from sqlserver to the ArcGIS server, that table (from sqlserver) is updating hourly, and the types are increasing.

I want to create dynamic layers according to type…( I am integrating those layers in Java script API)

So how should I do that? Any ideas/help will be highly appreciated..Thanks in advance : )

We have Desktop 10.1 and Server 10.1 only no ArcSDE.


texts:
I want to select raster layers by a GUI, so I designed a form by Qt designer that can add raster layers in 5 combo boxes. I want when a user select first raster layer in first combo box, and then user wants to select next layer in second combo box, name of first layer don't show in second combo box to select and this continue for other combo boxes.
my code is:

def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()

    # clear the combo box to list only current layers
    self.dlg.cBox_Exp_1.clear()
    self.dlg.cBox_Exp_2.clear()
    self.dlg.cBox_Exp_3.clear()
    self.dlg.cBox_Exp_4.clear()
    self.dlg.cBox_Sen_1.clear()
    self.dlg.cBox_Sen_2.clear()
    self.dlg.cBox_Capa_1.clear()
    self.dlg.cBox_Capa_2.clear()
    self.dlg.cBox_Capa_1.clear()
    self.dlg.linOutput.clear()

    # get the layers and add them to the combo box
    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.RasterLayer:
            self.dlg.cBox_Exp_1.addItem(layer.name(), layer)
            self.dlg.cBox_Exp_2.addItem(layer.name(), layer)
            self.dlg.cBox_Exp_3.addItem(layer.name(), layer)
            self.dlg.cBox_Exp_4.addItem(layer.name(), layer)
            self.dlg.cBox_Sen_1.addItem(layer.name(), layer)
            self.dlg.cBox_Sen_2.addItem(layer.name(), layer)
            self.dlg.cBox_Capa_1.addItem(layer.name(), layer)
            self.dlg.cBox_Capa_2.addItem(layer.name(), layer)
    # Run the dialog event loop
    result = self.dlg.exec_()


texts: How to save each band in a new row with it's name, in PostGIS?

texts:
I have used PostGIS's raster2pgsql to import my multi-band raster data into my PostGIS database. Here is the code that can do this operation:
raster2pgsql -I -C -e -Y -F /home/jj_user/input.tif output_table_name | psql -U my_db_username -d my_db_name -h localhost -p 5432

Here is it's SQL auto-generated codes in pgAdmin4:
-- Table: public.modis_5_band

-- DROP TABLE public.modis_5_band;

CREATE TABLE public.modis_5_band
(
    rid integer NOT NULL DEFAULT nextval('modis_5_band_rid_seq'::regclass),
    rast raster,
    filename text COLLATE pg_catalog."default",
    CONSTRAINT modis_5_band_pkey PRIMARY KEY (rid),
    CONSTRAINT enforce_srid_rast CHECK (st_srid(rast) = 0),
    CONSTRAINT enforce_scalex_rast CHECK (round(st_scalex(rast)::numeric, 10) = round(231.656358264, 10)),
    CONSTRAINT enforce_scaley_rast CHECK (round(st_scaley(rast)::numeric, 10) = round(- 231.656358264, 10)),
    CONSTRAINT enforce_width_rast CHECK (st_width(rast) = 41),
    CONSTRAINT enforce_height_rast CHECK (st_height(rast) = 41),
    CONSTRAINT enforce_same_alignment_rast CHECK (st_samealignment(rast, '0100000000DCCA0BE300F56C40DCCA0BE300F56CC0E17A146EA77552410358CBAA7D8A4F41000000000000000000000000000000000000000001000100'::raster)),
    CONSTRAINT enforce_num_bands_rast CHECK (st_numbands(rast) = 6),
    CONSTRAINT enforce_pixel_types_rast CHECK (_raster_constraint_pixel_types(rast) = '{8BUI,8BUI,8BUI,8BUI,8BUI,8BUI}'::text[]),
    CONSTRAINT enforce_nodata_values_rast CHECK (_raster_constraint_nodata_values(rast) = '{0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000}'::numeric[]),
    CONSTRAINT enforce_out_db_rast CHECK (_raster_constraint_out_db(rast) = '{f,f,f,f,f,f}'::boolean[]),
    CONSTRAINT enforce_max_extent_rast CHECK (st_envelope(rast) @ '01030000000100000005000000E17A146EA775524175E439B6F0774F41E17A146EA77552410358CBAA7D8A4F41A8345DE8ED7E52410358CBAA7D8A4F41A8345DE8ED7E524175E439B6F0774F41E17A146EA775524175E439B6F0774F41'::geometry) NOT VALID
)

TABLESPACE pg_default;

ALTER TABLE public.modis_5_band
    OWNER to postgres;
-- Index: modis_5_band_st_convexhull_idx

-- DROP INDEX public.modis_5_band_st_convexhull_idx;

CREATE INDEX modis_5_band_st_convexhull_idx
    ON public.modis_5_band USING gist
    (st_convexhull(rast))
    TABLESPACE pg_default;

This creates a table like below:

input.tif file has 5 bands. The above code saves all of the bands in one row.
What I want is to save each band in a new row with it's name. Each band's name is a string from a to e.
The desired table should look like this:


texts: Is the MakeFeatureLayer temporary location the same as "in_memory"?

texts:
From the ArcGIS Desktop Help:

  The layer that is created by the tool is temporary and will not
  persist after the session ends unless the layer is saved to disk or
  the map document is saved.

However, it does not state where this temporary location is.

Is it the same as the "in_memory" location that Esri suggests using while geoprocessing?


texts: How to get categorized list of all entities in a particular area of Google Maps/Open Street Maps?

texts:
I have been asked by my Professor to get the list of all kinds of entities for a particular area.

Here by entities I mean - Names and(or) Addresses categorically resolved to

residential area,
restaurants,
service stations or
any other possible category.

I am new to this field so I did some google search on how to get the list of such data for a particular area but was unable to find useful.

Is there any API for getting such data?

I don't want any code, I want people to point me in the right direction/give me some references.


texts: I am looking for a source of weather data for Europe. I use http://mesonet.agron.iastate.edu/cgi-bin/wms/goes/conus_vis.cgi? and http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r.cgi? for the US, but have been unable to find a source for Europe. I saw a similar question here: Listing available online WMS services (Weather, Land Data, Place Names)?, but it does not answer my question.


texts: Reprojecting coordinates using gdalwarp and OpenLayers?

texts:
I have an interpolated raster image with EPSG:27700 (British National Grid), after I warp it using gdalwarp into EPSG:3857 (WGS84 Web Mercator (Auxiliary Sphere)):

gdalwarp -s_srs EPSG:27700 -t_srs EPSG:3857 input.tif output.tif

which works fine (GDAL 1.10.1), (including a little bit of twisting (black cells) in the output one which is understandable as I guess there must be a little bit of distortion).

Then, if I create a Web Map Service (wms) in Geoserver using the output file as datastore, when I display the wms in a web page, by using OpenLayers and Google Maps as background, the image is not perfectly located, as it is misplaced about 90-100 meters LEFT where it should be in the British National Grid projection.

My question is whether gdalwarp does not the right transformation and/or OpenLayers does its best and does not adjust it, or it is definitely impossible to transform 27700 into 3857 accurately.


texts: What is the theory behind graduated symbols in QGIS?

texts: QGIS offers 5 methods of classifying vector map data (quantiles, equal interval, etc). What is the theoretical base of these methods? If there are certain methods suitable for certain types of data or so. Or if you know some webpages dealing with this. QGIS user guide just names them, doesn't explain.


texts: Extract Longitude and Latitude from Geotiff (ASTER 30m DEM)

texts:
my problem is apparently quite simple, but I cannot really find an easy answer: in QGIS I've created a raster layer by merging different GeoTiff images from the ASTER 30m DEM dataset.

The merging is ok: by moving the cursor I can see that the coordinates make sense and so the elevation.

What I'd like now is to extract the longitude and latitude informations from each pixel (generating a matrix the size of the layer with Lat values and one with Lon values), but I really don't know how to do it.

I also tried to import the image in MATLAB using the geotiffread function, but I cannot extract such information...

Thanks to all!


texts: more than one reactive events to subset markers in Leaflet in R

texts:
I am trying to create a shiny ap in R that displays markers on a map. I have a couple of categories that I would like to subset. Currently, I can subset but the map is only reactive on a single input (ethnicity). I would like the map to be reactive no matter the input. I tried doing the following

  df <- eventReactive(input$Ethnicity | input$IQ, {
     daf <- subset(data_F, (Ethnicity == input$Ethnicity & IQ > as.numeric(input$IQ)))
     subset(daf, select = c("longitude", "latitude"))
  })

This unfortunately did not work. So currently I would like an Or event statement. Either a user can manipulate the slider or the pull down menu and update the markers on the map.

Currently I am failing drastically.

Here is my code as it stands.

ui <- fluidPage(
                selectInput(inputId = "Ethnicity",
                            label = "Ethnicity",
                            choices = c("Green","Purple")),

            sliderInput (inputId = "IQ",
                    label = "IQ", step = 1,
                    value = 0, min = 0, max = 200)
                    ,
  leafletOutput(outputId = "map")
)

server <- function(input, output, session){

  df <- eventReactive(input$Ethnicity, {
     daf <- subset(data_F, (Ethnicity == input$Ethnicity & IQ > as.numeric(input$IQ)))
     subset(daf, select = c("longitude", "latitude"))
  })

  output$map <- renderLeaflet({
   leaflet () %>%
    addTiles() %>%
    addMarkers( lng = df()$longitude, lat = df()$latitude,
        clusterOptions = markerClusterOptions())
  })
}

shinyApp(ui = ui, server = server)


texts: ArcMap 10.4.1.5686. What I am asking may be impossible, but it sure would be sassy. I went to a site and GPS marked areas where people had been creating ATV trails. I want to link photos of the damage to the points in such a way that if you mouse over the point, a photo pops up. I have tried the adding photos to attribute table thing, but when I click identify, info just pops up, and then I have to click again to see the photo. I'd like to streamline this if at all possible.


texts: Map project a raster having separate latitude and longitude raster bands

texts:
I have a satellite image where one of the channels is a raster having the longitude values, another channel a raster containing the latitude values.

What would be the best way for map projecting such files? Is there some standard routines or tools (I use Python / gdal preferrably)?

My only thought so far is to loop through the pixels, convert the lat/long value of a given pixel and drop it into the closest matching pixel of a raster already in the projection I want to have (where it seems not that straightforward to actually find the closest matching value) -- but somehow I wonder if there is a better and quicker way?


texts:
I need to build a solution that allows me to track cellphones/tablets and display on a visual map.

I think there are two parts, one is understanding the general architecture. The other is more detail. I have seen a few websites that do what i want to do.For example,the Android Device Manager lets you see where your devices are on a map. I'd envisage the tablets to be Android based so even if there is a SAAS solution that I could plug into my web page would be perfect

Basically I was thinking of installing a generic app(possibly android) on a tablet/cellphone then capturing the gps data somehow through a web interface and displaying the tablets onto a single map on my webpage. A SAAS solution would work well or alternatively a javascript package to read the captured GPS co ordinates and display on a map.

Any help is appreciated.
Regards


texts:
I am having a Merging issue with the GDAL python library. I am using gdal_merge.pylink
My issue is when I merge two rasters NoData value of one raster overlaps the other data resulting in not a proper raster image.
For example, in the raster images below (left images) if I merge these 4 images labeled as 1,2,3,4 I should get the merged raster image as the right-sided image.

But, the merged image I am getting looks like this (image below):
(Note: Image labeled 1,2 and 3 are merged properly. It's only when we merged them with image labeled 4, it's not properly merging.)

If anyone knows, how to resolve this issue then it'll be a great help to me. Thanks in advance.


texts: How to output True/False if attribute is a number in QGIS field calculator?

texts:
I have a points shapefile which has two fields that I want to use to determine the size of each point. One field, "StepYear", refers to the year of step change in a particular time series, and it contains both numbers (e.g. 1996), and non-numericals (such as NA or - ). The non-numericals indicate a number of different cases in which the test did not proceed.

If the test did not proceed, then I want the size of the point to be equal to 2. Otherwise, I want the size to be determined by the statistical significance of the test. I want code similar to the following:

CASE
  WHEN  NOT ISNUMBER("StepYear") THEN 1
  WHEN  "StepSig" <= 0.01 THEN 8
  WHEN  "StepSig"  <= 0.05 THEN 6
  WHEN  "StepSig"  <= 0.1 THEN 4
  ELSE 2
END

However, I haven't been able to find any function that performs the "ISNUMBER" check. Is there a function in the QGIS field calculator which outputs TRUE/FALSE for a numerical type check?


texts: Getting elevation for multiple lat long coordinates in Python

texts:
I have seen the answer to Elevation info for multiple lat long coordinates and I would like to be able to do this with a python script instead. I cannot use Googles API because I have 8000 points. I have them as part of a .csv file which I import into pandas as a dataframe.

something like:

elevation_list = []

def elevation_function(df):
    for lat,long in zip(df.lat,df.lon):
          elevation_i = some_call_to_file(lat,long)
          elevation_list.append(elevation_i)

My coordinates span the western half of the US so Im not sure the function could be so simple. I have downloaded all SRTM90M resolution tiffs from this amazing person


texts: Adding points randomly inside a polygon from an attribute table to a different polygon layer

texts:
I have two layers:
1. One with polygons geometry
2. CSV table with more than 600k rows

They have an atribute ID in common, so, I want to distribute randomly each row, inside the polygon as a point. I'm using QGIS 3.4


texts: R lidR does not read the CRS of my las file

texts:
I have a LAS 1.4 file that I can read in lidR but the CRS is not read:
las = readLAS("LD31201353.las")
las
#> class        : LAS (v1.4 format 6)
#> memory       : 172.5 Mb
#> extent       : 3121464, 3123000, 1353000, 1356000 (xmin, xmax, ymin, ymax)
#> coord. ref.  : NA
#> area         : 4.5 kunits²
#> points       : 1.97 million points
#> density      : 0.44 points/units²

If I look at the header the VLR is empty
las@header
#> File signature:           LASF
#> File source ID:           0
#> Global encoding:
#>  - GPS Time Type: Standard GPS Time
#>  - Synthetic Return Numbers: no
#>  - Well Know Text: CRS is WKT
#>  - Aggregate Model: false
#> Project ID - GUID:        00000001-0012-07e3-0000-000000000000
#> Version:                  1.4
#> System identifier:        Merrick
#> Generating software:      MARS
#> File creation d/y:        188/2018
#> header size:              375
#> Offset to point data:     377
#> Num. var. length record:  0
#> Point data format:        6
#> Point data record length: 30
#> Num. of point records:    1965664
#> Num. of points by return: 1634171 285040 42901 3469 82 1 0 0 0 0 0 0 0 0 0
#> Scale factor X Y Z:       0.001 0.001 0.001
#> Offset X Y Z:             3122231 1354500 11639
#> min X Y Z:                3121464 1353000 11008.84
#> max X Y Z:                3123000 1356000 11782.48
#> Variable length records:  void

But the file actually contains a CRS because lasinfo from LAStools can read it
extended variable length header record 1 of 2:
 reserved             0
 user ID              'LASF_Projection'
 record ID            2112
 length after header  1213
 description          'OGC COORDINATE SYSTEM WKT'
   OGC COORDINATE SYSTEM WKT:
   COMPD_CS["NAD83(2011) / Colorado Central (ftUS) + NAVD88 height (ftUS) - US Geoid Model of 2012 B",PROJCS["NAD83(2011) / Colorado Central (ftUS)",PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["central_meridian",-105.5],PARAMETER["latitude_of_origin",37.833333333333],PARAMETER["standard_parallel_1",39.75],PARAMETER["standard_parallel_2",38.45],PARAMETER["false_easting",3000000],PARAMETER["false_northing",1000000],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6428"]],GEOGCS["NAD83(2011)",DATUM["NAD83 (National Spatial Reference System 2011)",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["US survey foot",0.3048006096012192, AUTHORITY["EPSG","9003"]],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","6318"]],VERT_CS["NAVD88 height (ftUS) - US Geoid Model of 2012 B",VERT_DATUM["North American Vertical Datum 1988",AUTHORITY["EPSG","5103"]],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","6360"]],PARAMETER[Vertical_Shift,0.0],PARAMETER[Direction,1.0]]

What do I need to do in R to expose the CRS (besides calling lastools from R...)?


texts:
To give context, I was using the Intersect Tool to find the number parcels that are located in multiple Sections of an STR.  I tried both the Geoprocessing Tool, and calling the function in ArcPy.  I noticed I got different results when my input was the same order for both methods, Parcel Layer first, then Section layer. However when I did  Parcel Layer first, Section Layer second for the Geoprocessing Intersect Tool, then Section first, Parcel Second with the ArcPy Function, I obtained the same results.

Does order matter with the Intersect Tool?

I figure either it does, or I was just doing something wrong in my ArcPy script that was giving inconsistent results.


texts: NOAA API data for past monthly mean Temperature and precipitation

texts:
I am trying to use the NOAA API to get data on monthly mean temperature and precipitation for the past 10 years. I wrote the following python code (the variable "token" in the code is the API token that I got:

import requests
url = "http://www.ncdc.noaa.gov/cdo-web/api/v2/data?datasetid=GHCNDMS&locationid=ZIP:28801&startdate=2000-01-01&enddate=2010-01-01"
response = requests.get(url, headers = token)
data = response.json()

I am not able to get data for for the specified time period. I am not sure where I am going wrong.


texts: Overpass Turbo: query for administrative level results in empty dataset

texts:
If I look up a city on openstreetmap, I can see it has an administrative level of 10. For example: https://www.openstreetmap.org/relation/162798

Subsequently, I run a query in Overpass Turbo to extract the aministrative boundary of this city:

[out:json][timeout:60];
{{geocodeArea:'s-Hertogenbosch}}->.searchArea;
    relation["admin_level"="10"](area.searchArea);
);

out body;
>;
out skel qt;

{{style:
  area
  { color:gray; fill-color:DarkGray; }
}}

However, this does not work (returns an empty dataset). If I set the admin_level to 8, it returns a much bigger area compared to what openstreetmap shows. For other cities, the admin_level 10 does work, so why not for 's-Hertogenbosch? Am I missing something here?


texts: How to display survey calls using the traverse tool?

texts: I use the Traverse Tool to run out legal descriptions on deeds and recorded subdivision plats. I'm looking for a way to display the survey calls once I'm done. Is there a way to do that? I'm using ArcMap 10.0.


texts:
I'm having issues with my portal enterprise 10.6 setup.

I have a feature service which was created using Survey123. It all worked fine but one day all attachments have gone missing and if I checked the logs I would get the following:

I exported the Feature service as an FGDB and it all looks fine inside a GDB.

I would just remake the service and move the data but I had this issue in the past and I would like to find a solution rather than a workaround.


texts: osm2pgsql Error: operator class "gist_geometry_ops" does not exist for access method "gist"

texts:
When using osm2pgsql to import massachusetts-latest.osm.bz2 into PostGIS 2.0.1 database osm (running on PostgreSQL 9.1 on an Ubuntu 12.04), I get the error failed: ERROR:  operator class "gist_geometry_ops" does not exist for access method "gist".

I have ran the enabler .sql scripts including legacy.sql. Truncating the tables and reimporting does not help. Any idea what is causing this error?

Command

osm2pgsql massachusetts-latest.osm.osm.bz2 -d osm -U postgres -P 5432 --hstore -W

Output:

Reading in file: massachusetts-latest.osm.bz2
Processing: Node(23054k) Way(2370k) Relation(5058)  parse time: 300s

Node stats: total(23054254), max(2147478725)
Way stats: total(2370244), max(215072011)
Relation stats: total(5058), max(2855687)

Writing way (2370k)

Writing relation (5058)
Committing transaction for planet_osm_point
Sorting data and creating indexes for planet_osm_point
Committing transaction for planet_osm_line
Committing transaction for planet_osm_roads
Sorting data and creating indexes for planet_osm_line
Committing transaction for planet_osm_polygon
Sorting data and creating indexes for planet_osm_roads
Sorting data and creating indexes for planet_osm_polygon
CREATE INDEX planet_osm_roads_index ON planet_osm_roads USING GIST (way GIST_GEOMETRY_OPS) TABLESPACE pg_default;
 failed: ERROR:  operator class "gist_geometry_ops" does not exist for access method "gist"

Error occurred, cleaning up


texts: Up-sampling (increasing resolution) raster image using GDAL?

texts:
I have a single-band geo-referenced raster image (a DEM) and my goal is to increase the number of pixels in each dimension (x and y) by 3.
So given an initial resolution of 1200 * 1200 pixels, I want to upsample it to 3600 * 3600 pixels, keeping the geographic referencing.

This procedure is a part of a large batch processing so simple gui workflow is not an option and I'd like to do it programmatically, using python gdal.

What is the right way to do it and is there any interpolation required to do it?


texts: Maximum number of point features in an OpenLayers vector layer

texts:
In your experience, how many point features can be added to an OpenLayers vector layer (new OpenLayers.Layer.Vector("Point Layer")) before it goes unusably slow?

My use case is to display points from a database table. The user can decide which time frame to visualize. Therefore the result can be from very few to potentially 100,000s of points. I'd like to introduce a reasonable limit and warn the user if his query would return more features.


texts: Return latitude and longitude from geometry point type

texts:
For storing latitude and longitude in the database, I use the following query:

WITH cte_posttype AS (
 SELECT posttype_id, code AS posttypeCode FROM "post"."posttype" WHERE code = $6;
)
INSERT INTO "post"."post"(post_id, user_id, fb_username, place_name, geom, posttype_id, status_text, status_img)
VALUES ($1, $2, $3,
   ST_SetSRID(ST_MakePoint($5, $4), 4326),
   SELECT posttype_id FROM cte_posttype,
   $7, $8)
RETURNING post_id, user_id, fb_name, place_name, geom, cte_posttype.posttypeCode, status_text, status_img, updated_at, created_at;

The relevant part is:

ST_SetSRID(ST_MakePoint($5, $4), 4326)

Now how do I revert the geometry type back to get the original lat and lng to return it?
The same question applies on making a SELECT: I want to retrieve the lat and lng value to return it to the client as such.


texts: Help required to convert georeferenced (and rectified) TIFF orthophotos to georeferenced .JPG format

texts:
I am not a GIS specialist (so patience required!) :) We work with georeferenced, rectified 1:10 000 orthophotos. I need help with the conversion of a georeferenced tiff to a georeferenced jpg.

Our recently acquired JUNO (Trimble) GPS's work very well with an eventual conversion of jpg to ecw format - but I first need help in converting my tiff files to jpg.

Any ideas please? Free/shareware or products for sale that you know of?
thanks!
Monique


texts: Error-'ERDAS Imagine could not find sensor information(NDVI calculation)'

texts: I want to calculate NDVI from a Landsat 8  image using ERDAS imagine 2020. However the software is giving me an error message that reads 'Warning: Could not find sensor information...' for the image i have used(see screenshot).


texts: Use multiple filter for Tiled layer in GeoServer

texts:
I am using GeoServer & want to apply multiple filters for 'tiled layer'. Is it possible?

I checked this documentation & I have researched on web regarding this but unable to find exact documentation


texts: Google Earth Engine: Dropdown of items in ImageCollection

texts:
I want to create a select menu which lists all datasets in an ImageCollection. When the user selects an item, the referenced dataset should be displayed in the map. Here's a static example that works and illustrates my approach so far...

Step 1 - Acquire L8 scenes for region of interest and create a list of what should be the label-value pairs for the select menu:

var roi = ee.Geometry.Rectangle(113, 53.5, 114, 54);
Map.centerObject(roi);
Map.addLayer(roi);

var landsatCollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
    .filterBounds(roi);

var sat_scene_dict = landsatCollection.map(function(sat_scene) {
  var dict = {
    value: sat_scene.id(),
    label: sat_scene.get('DATE_ACQUIRED'),
  };
  return ee.Feature(sat_scene.geometry(), dict);
});

Step 2 - Initialize select items and populate with values

var imageSelect = ui.Select({
  items: [ // Replace this part with labels/value from the list of scenes
    {label: '2013-06-16', value: 'LC08_128022_20130616'},
    {label: '2015-05-21', value: 'LC08_128023_20150521'}
  ],
  onChange: function(value) {
    var selected_scene =
        landsatCollection.filter(ee.Filter.eq('system:index', value));
    Map.addLayer(selected_scene, {bands: ['B3', 'B2', 'B1']}, 'Selected scene');
  }
});

So, the whole problem is about how to populate the 'items' of ui.Select with label-value pairs from the list of scenes. As soon as I put any form of sat_scene_dict to items I get errors in in the form of The argument provided to ui.List.reset() must be an array.. I have tried to convert the sat_scene_dict to a list of dictionaries as in the static version but keep getting errors that map must return features or images.


texts:
I have got two geometries one is polygons (postcode boundaries) and second is linestring (roads). I want to intersect them together and want roads that are intersecting the polygons. While doing this I got that many sections of roads are  still connected to each other (dissolved). So I want to separate them (explode) as they are in geometry. I did search for the solution but haven't found anything.

This is the case below they all are a one feature in QGIS. I want them as single entity.


texts: PostGIS geography query returns a string value

texts:
I have a strange issue. The lonlat column on my app works well on the development server –– its output is in the form of POINT(X Y). But when I move the data to the production server, the output is strange!

ActionView::Template::Error (undefined method `lon' for "0101000020E6100000541B9C887E7A52C02920ED7F80614440":String):

The lonlat value, which is encoded with SRID: 4326, is being read as a string. I am almost certain that there was a corruption in the data during migrating it from development to production because this was not a problem before the migration.

Does anyone know what about the database schema or column may cause this issue?


texts: How do I create a basic image layer as a base layer in OpenLayers?

texts:
I'm trying to create a map using OpenLayers which uses a png image as a base layer but I'm having a bit of trouble. I'm not too familiar with OpenLayers as I have only come across it today and thought that it could prove to be a powerful tool.

So far I have just been using the code from the tutorial and changed the base layer from a WMS format to an Image format but it still does not show the map on the webpage.

<!DOCTYPE html>

    <title>OpenLayers Tutorial - Basic Map Setup</title>

OpenLayers Tutorial - Basic Map Setup

It's probably something really silly that needs changing but I've looked at documentation and I can't seem to find the issue.

Any help would be appreciated!


texts:
I am trying to generate the maritime distance between seaports across the world. I am using the cost distance and cost path tool from ArcGIS to extract the shipping lanes from a raster file which contains the density of ships.

Here is the problem that I have. When I tried finding the path connecting Shanghai and San Francisco it provides the result in the figure below. The red line is the path found, the green circle is Shanghai and the green square is San Francisco. As can be seen, it takes all the trouble to go through the Atlantic and Indian Ocean, instead of just cross the Pacific.

I think the problem is that the two sides of the Pacific are not connected on the map. Is there a way that I can tell the tool that the left hand side and right hand side of the map are connected? The projection that I am using is WGS_1984_Mollweide.


texts: Earth Engine: sampleRectangle 262144 pixel limit

texts:
I am attempting to use the ee.Image.sampleRectangle() to isolate a selected region. However, I obtain the error EEException: Image.sampleRectangle: Too many pixels in sample; must be <= 262144. Got 372096.
I realize that there is a limit of 262144 pixels for ee.Image.sampleRectangle(). I have heard that I can use my google drive to overcome this limit, but I am not sure on how to go about that.
Edit:
The exact code snippet that I used was
rgb_img = areaofearth.sampleRectangle(defaultValue=0)
rgb_img = np.array(rgb_img.getInfo()['properties']['b1']) #This line gives the error
rgb_img = rgb_img.reshape(len(rgb_img), -1)

I am essentially trying to convert the ee.Image titled areaofearth to an image via numpy array. The above code works great for small areas.


texts: Transforming regular lat/lon to rotated lat/lon?

texts:
I am currently working with climate data projected in a rotated lat/lon grid and came across this excellent 'cookbook' by @mkadunc on Manually transforming rotated lat/lon to regular lat/lon?

I’ve implemented the solution for transforming a rotated grid to a regular grid; however, I’m also interested in doing the opposite. My first thought was to just change the operational sign of ϑ and φ, so that the transformation will go ‘the other way around the clock’, but this doesn’t seem to work leaving me a bit stuck here. I’ve spend some hours now trying to read up on my trigonometry, but still haven’t been able to ‘reverse’ the function.


texts: A doubt on the algorithm Cross-Track Distance in Java Code

texts:
I was researching a way to calculate the shortest distance between a coordinated and a straight when I met the Cross-Track Distance. To find the distance, I had to implement the following code:

final int EarthRadius = 6371; //in kilometers
double lon1, lat1, lon2, lat2, lon3, lat3;

double brng13 = bearing(lat1, lon1, lat2, lon2); // in degrees
double brng12 = bearing(lat1, lon1, lat3, lon3); //  in degrees
double dist13 = harvesine(lat1, lon1, lat3, lon3); // in kilometers

double crosstrack = Math.asin(Math.sin(dist13/EarthRadius)*Math.sin(Math.toRadians(brng13 - brng12))) * EarthRadius;

As I'm implementing the code in Java, functions like sine and cosine receive as argument an angle in radians. Therefore, the line of the cross track should be like this:

double crosstrack = Math.asin(Math.sin(Math.toRadians(dist13/EarthRadius))*Math.sin(Math.toRadians(brng13 - brng12))) * EarthRadius;

The problem is that during my research, I found several codes in Java without converting "dist13/EarthRadius" to radians.

What is the correct way? Convert to radians or not?

Edited - Added missing java functions:

public static double haversine(double lat1, double lon1, double lat2, double lon2) {
        double dLat = Math.toRadians(lat2 - lat1);
        double dLon = Math.toRadians(lon2 - lon1);
        lat1 = Math.toRadians(lat1);
        lat2 = Math.toRadians(lat2);

        double a = Math.pow(Math.sin(dLat / 2),2) + Math.pow(Math.sin(dLon / 2),2) * Math.cos(lat1) * Math.cos(lat2);
        double c = 2 * Math.asin(Math.sqrt(a));
        return EarthRadius  * c;
    }
public static double bearing(double lat1, double lon1, double lat2, double lon2){
  double longitude1 = lon1;
  double longitude2 = lon2;
  double latitude1 = Math.toRadians(lat1);
  double latitude2 = Math.toRadians(lat2);
  double longDiff= Math.toRadians(longitude2-longitude1);
  double y= Math.sin(longDiff)*Math.cos(latitude2);
  double x=Math.cos(latitude1)*Math.sin(latitude2)-Math.sin(latitude1)*Math.cos(latitude2)*Math.cos(longDiff);

  return (Math.toDegrees(Math.atan2(y, x))+360)%360;
}


texts: Combining multiple shapefiles into one shapefile using ArcMap?

texts:
I have two .shp containing stream segments from two different NHD HUC files (see photo). I would like to make these into one shapefile and combine the attribute tables. I have tried to use the merge tool in ArcMap but the attribute tables do not match like they were before being merged. Am I missing a step?

I am using ArcMap 10.5 and I would like the attributes of one shapefile within the same table as the attributes of the other shapefile in one table as one shapefile.


texts:
I’m new with Arc 10, and I'm struggling alot.
When using Geoprocessing I’m getting the same error or message again and again.
"ERROR 999998: Unexpected Error"

I have checked all the settings I can think off.
Any help or recommendation will be greatly appreciated.


texts: FME SDE Writer changing part sequence of Multi-part polygon on Insert?

texts:
I am creating a FME script that will accept to databases (MDB file and SDE Database) and will update the SDE database if there are any changes to the MDB file. I have a CRCCalculator for both SDE and MDB file to determine if there are any changes and a UniqueIdentifier is present on both.

Its currently working but I'm noticing that when FME writes the features to the SDE database (e.g. FMEMultiArea), in some cases it changes the sequence of the parts therefore changing the area value. This will in turn change the CRC value calculated when i run the FME script again so everytime it thinks there is an update.

I'm wondering if i need to add something to prevent the writer from changing the sequence.

FME Workflow:

Sample MDB Feature:

Sample SDE Feature:

I am using SQL Server 2012 to store the data.


texts: Looking for atmospheric dispersion model with 3d GIS application

texts: is there any atmospheric dispersion model with 3d GIS application that can use DEM and atmospheric data as input ?


texts:
Using mapserver 7.0, is there a way to suppress the error messages sent to the user, such as :-

<?xml version='1.0' encoding="UTF-8" standalone="no" ?>
  <!DOCTYPE ServiceExceptionReport SYSTEM "http://schemas.opengis.net/wms/1.1.1/exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1">
<ServiceException>
  msDrawMap(): Image handling error. Failed to draw layer named &#39;FOO&#39;.
  msPostGISLayerOpen(): Query error. Database connection failed. Check server logs for more details.Is the database running? Is it allowing connections? Does the specified user exist? Is the password valid? Is the database on the standard port?
 </ServiceException>
</ServiceExceptionReport>

Running in the apache environment.

I have tried setting the ERROR and EMPTY urls in the WEB header but they do not seem to work.


texts: Shortest Path between large area -- RAM issue

texts:
I built a shortest path algorithm that finds and stores the distance of the optimal path between neighboring grid cells using a friction surface as an input. It works beautifully but I am facing RAM memory issues to do so for a very large area, which is my ultimate goal. My object of interest is tau, the bilateral distances' matrix below.

library(raster)
library(rgdal)
library(rgeos)
library(gdistance)
library(geosphere)

load('grid2.rdata') # small grid of locations whose bilateral path I need to calculate
# load('grid.rdata') # very large grid which I would like to use
grid <- area
# friction surface:
friction <- raster('test.tif')
# standardizing projections:
grid <- spTransform(grid, friction@crs)
# cropping friction into aoi area:
friction <- crop(friction,grid)

# making it a transition matrix:
friction <- transition(friction, function(x) 1/mean(x), 8)

# Now identifying which cells are neighbors of the others:
tau <- as.matrix(gTouches(grid, byid = T))
rownames(tau) <- colnames(tau) <- 1:nrow(grid)
# tau is symmetric -- replace upper triangle
# part with 0 (cuts iterations by half!):
tau[upper.tri(tau)] <- 0
ll <- sum(tau) # number of total iterations

# Calculating routes between neighboring cells:
ii<-0 # counter
cat('2. Calculating shortest route between cells:\n')
for (i in 1:nrow(grid))  {
  for (j in which(as.numeric(tau[i,])==1)) {
    ii <- ii+1
    pp <- SpatialPoints(coords = matrix(c(grid@data[i,long],grid@data[i,lat], grid@data[j,long], grid@data[j,lat]),2,2,byrow=T), proj4string = grid@proj4string)
    d <- shortestPath(friction,
                      pp[1,],
                      pp[2,],
                      output = "SpatialLines"
    )
    # Desired outputs:
    if (ii==1) {routes <- d} else {routes <- rbind(routes, d)}
    if (ii==1) {locs <- pp} else {locs <- rbind(locs, pp)}
    # storing the distance (in km):
    tau[i,j] <- lengthLine(d)/1000

  }

}

It works pretty well and if I want to visualize the output with

plot(grid)
plot(locs,pch=T,add=T)
plot(routes,add=T)

I obtain the image attached. How could I do the same for a very large area (grid.rdata at the beginning of the code) in a more memory-efficient way?

link to the data used: friction surface, grid, grid2.